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ABSTRACT 

Redundant  Spacing  Calibration  (RSC)  techniques  employ  redundancy  in  the  baselines 
of  a  telescope  array  to  eliminate  the  contribution  of  atmospheric  turbulence  in  the  in¬ 
terferometric  observables.  Whereas  conventional  techniques  for  this  elimination  require 
the  enforcement  of  prior  constraints  on  the  underlying  image,  RSC  algorithms  can  be, 
in  principle,  mathematically  well-posed  and  hence  require  no  such  prior  knowledge. 
Traditionally  these  algorithms  have  been  applied  directly  to  the  fringe  measurements. 
However,  in  scenarios  of  low  photon  flux,  such  as  those  arising  in  the  observation  of 
dim  objects  in  space,  single-exposure  fringe  measurements  are  not  reliable  observables 
in  general.  Instead  one  must  rely  on  time-averaged,  atmosphere-invariant  quantities 
such  as  the  bispectrum.  In  this  paper,  we  develop  a  novel  algorithm  for  redundant  ar¬ 
rays  which  provides  robust  image  reconstruction  using  integrable  atmosphere-invariant 
observables.  Our  algorithm  utilizes  standard  linear  estimation  methods,  as  well  as 
techniques  from  lattice  theory,  to  reliably  estimate  the  Fourier  phase.  Moreover,  we 
provide  theoretical  and  empirical  evidence  that  generalizing  the  classical  bispectrum 
to  higher-order  atmosphere-invariant  observables,  which  we  call  n-spectra,  can  offer 
significant  performance  gains.  Our  selection  of  an  independent  and  higli-SNR  set  of 
n-spectra  leverages  the  notion  of  the  minimum  cycle  basis  from  graph  theory.  We  an¬ 
alyze  the  expected  shot-noise-limited  performance  of  our  algorithm  for  both  pairwise 
and  Fizeau  interferometric  architectures,  and  corroborate  this  analysis  with  simulation 
results  showing  performance  near  the  Cramer-Rao  bound.  Lastly,  we  apply  techniques 
from  the  field  of  compressed  sensing  to  perform  image  reconstruction  from  the  esti¬ 
mated  complex  visibilities. 

Key  words:  techniques:  interferometric  -  techniques:  image  processing  -  atmospheric 
effects  -  methods:  analytical  -  methods:  numerical 


1  INTRODUCTION 


This  measurement  is  typically  obtained  via  compu¬ 
tation  in  radio  interferometers,  since  at  radio  wave¬ 
lengths  the  fields  can  be  sampled  and  correlated  dig¬ 
itally.  Optical  interferometers,  on  the  other  hand, 
obtain  the  correlation  indirectly,  by  measuring  the 
intensity  of  the  interference  fringe,  i.e.  superposi¬ 
tion  of  the  fields  incident  at  each  aperture.  As  es¬ 
tablished  by  the  van-Cittert  Zernike  Theorem  (van 
Cittert  1934)  (Zernike  1938),  each  correlation  pro¬ 
vides  a  sample  of  the  scene’s  2D  Fourier  Transform. 
Over  the  past  few  decades,  a  diverse  set  of  algo¬ 
rithms  has  been  developed  to  reconstruct  an  image 
of  the  scene  from  these  samples. 


Interferometers  can  produce  high-resolution  im¬ 
agery  of  a  distant  object  by  combining  signals  col¬ 
lected  by  an  array  of  apertures.  Assuming  the  spa¬ 
tial  extent  of  the  array  is  sufficiently  large,  they 
can  in  principle  achieve  the  resolution  of  a  sin¬ 
gle  aperture  of  prohibitively  large  size.  The  critical 
measurement  enabling  image  reconstruction  is  the 
time-averaged,  pairwise  cross-correlation  between 
the  fields  observed  at  each  aperture  in  the  array. 
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A  principal  challenge  in  interferometric  image  recon¬ 
struction  is  variation  in  the  complex  gains  among  the  multi¬ 
ple  apertures  of  the  interferometer.  This  variation  can  arise 
from  differences  among  the  analog  components  of  the  an¬ 
tenna  elements  (e.g.  cable  length  differences)  in  the  array. 
Moreover,  atmospheric  turbulence  alters  the  effective  path 
length  (or  optical  piston)  from  the  target  to  each  telescope 
by  a  non-uniform  and  time- varying  amount.  For  simplic¬ 
ity,  we  will  henceforth  refer  to  such  aperture-specific  phase 
variation  as  optical  path  difference  (OPD),  regardless  of  its 
source.  As  a  result  of  OPD,  the  Fourier  component  mea¬ 
sured  by  apertures  k  and  j  is  given  by  y%;  =  \ykj\e‘^kj  ■  in 
which  =  (0kj  +  S4>kj)  m°d  2tt ,  and  is  the  true  Fourier 

phase  at  spatial  frequency  ,  and  S (p^j  represents  the  OPD 
observed  between  apertures  k  and  j.  In  a  non-redundant  ar¬ 
ray  of  N  apertures,  there  are  such  measurements  cor¬ 
responding  to  each  pair  of  apertures.  Collecting  such  mea¬ 
surements  together,  we  obtain  a  system  of  equations: 

P  =  0+A<p  (1) 

where  A  is  an  x  N  matrix  whose  rows  contain  the  piston 
differences  involved  in  each  measurement. 

Since  there  are  unknown  Fourier  phases  and  N  -  1 
unknown  OPDs,  inference  of  the  Fourier  phases  from  the 

phase  measurements  is  inherently  ill-posed  for  such  a 
non-redundant  array.  The  traditional  means  of  mitigating 
this  issue  has  been  to  impose  additional  prior  constraints  on 
the  image  to  compensate  for  this  degeneracy.  For  example, 
many  popular  methods  assume  and  enforce  smooth¬ 
ness  of  the  image  by  bounding  its  total  variation 1 
(Rudin  et  al.  1992).  Other  effective  approaches  fa¬ 
vor  a  pixel  distribution  of  high  entropy  (Skilling  & 
Bryan  1984). 

Reconstruction  in  optical  interferometry  amounts  to 
finding  an  image  within  prior  constraints  that  also  agrees 
with  the  observed  data.  To  achieve  this  objective,  a  myriad 
of  so-called  self-calibration  algorithms  were  developed  in  ra¬ 
dio  astronomy  (see,  e.g.,  Pearson  &  Readhead  (1984)  and 
Besnerais  et  al.  (2008))  which  alternate  between  estimation 
of  the  Fourier  phases  and  estimation  of  the  image  coefficients 
subject  to  the  imposed  constraints.  Such  algorithms  have 
shown  to  be  successful  in  high-flux  scenarios  in  which  the 
interferometric  measurements  have  a  high  signal-to-noise  ra¬ 
tio  (SNR).  In  low-flux  scenarios  on  the  other  hand,  one  must 
integrate  observables  over  a  long  time  period  in  order  to 
build  sufficient  SNR  for  reliable  image  reconstruction.  Since 
integration  beyond  the  coherence  time  of  the  atmosphere 
would  result  in  fringe  blurring,  one  is  then  led  to  the  forma¬ 
tion  and  subsequent  integration  of  atmosphere-invariant  ob¬ 
servables  from  each  atmosphere-coherence  time  (or  frame). 
The  classic  of  these  observables  is  the  triple  product  of  the 
Fourier  components  along  the  sides  of  a  baseline  triangle 
(e.g.  b  12, 1>23>  and  l>3i).  Note  that  the  OPD  cancels  in  these 
products  and  hence,  like  the  Fourier  magnitudes,  these  so- 
called  bispectra 2  are  OPD-invariant  observables.  However, 

1  The  total  variation  of  an  image  is  defined  as  the  integral  of  the 
absolute  value  of  the  image  gradient. 

“  It  is  customary  in  the  literature  to  use  bispectrum  to 
refer  to  the  entire  function  multiplying  Fourier  samples 


for  a  non-redundant  array  with  ^ )  distinct  baselines,  re¬ 
covery  of  the  Fourier  phases  from  the  bispectra  phases  (i.e. 
the  closure  phases)  remains  ill-posed  since  there  are  only 
independent  bispectra  (Readhead  et  al.  1988).  Suc¬ 
cessful  bispectra-based  image  reconstruction  remains  feasi¬ 
ble  in  spite  of  this  ill-posedness  (see  e.g.  Thiebaut  (2013), 
Besnerais  et  al.  (2008)),  but  again  prior  constraints  (e.g. 
on  the  image  support)  must  be  enforced  to  regularize  the 
reconstruction.  It  is  well-known  that  such  prior-regularized 
approaches  become  unreliable  when  the  scene  under  obser¬ 
vation  does  not  match  the  prior  assumptions  being  enforced. 

As  interferometric  systems  are  tasked  with  the  observa¬ 
tion  of  increasingly-complex  scenes,  approaches  which  mini¬ 
mize  reliance  on  specialized  prior  assumptions  are  expected 
to  gain  appeal.  An  alternative  and  intrinsically  well-posed 
approach  to  prior-regularized  phase  recovery  is  to  use  base¬ 
line  redundancy  to  explicitly  solve  for  OPD  variation;  an 
array  with  baseline  redundancy  contains  repeated  instances 
of  the  same  baseline  involving  distinct  aperture  pairs.  Since 
Fourier  phases  can  be  assumed  to  be  equal  for  all  repeated 
baselines,  an  observed  difference  amongst  their  correspond¬ 
ing  measurements  exposes  the  contribution  of  the  OPD.  This 
idea  of  using  redundant  arrays  to  calibrate  out  OPD  varia¬ 
tion,  known  as  Redundant  Spacing  Calibration  ( RSC ),  was 
developed  in  works  such  as  those  by  Hamaker  et  al.  (1977), 
Arnot  et  al.  (1985)  and  Greenaway  (1990).  In  recent  years, 
innovation  in  optical  technology  has  engendered  a  revival  of 
interest  in  the  RSC  technique.  The  simultaneous  (or  Fizeau- 
style)  measurement  of  fringes  on  a  common  focal  plane  has 
long  been  a  popular  method  of  acquiring  many  baseline  mea¬ 
surements  in  an  economical  manner.  However,  the  Fizeau 
method  had  been  incompatible  with  RSC  techniques  since 
the  fringes  formed  by  each  set  of  redundant  baselines  would 
alias  on  the  focal  plane.  An  elegant  solution  to  this  problem 
was  proposed  by  Perrin  et  al.  (2006).  This  work  developed 
the  idea  of  segmenting  the  entrance  pupil  of  a  single  tele¬ 
scope  into  an  RSC  arrangement  of  sub-pupils  from  which 
the  light  was  then  coupled  via  single-mode  fiber  to  a  non- 
redundant  exit  pupil,  thereby  permitting  unambiguous  and 
simultaneous  fringe  detection  for  an  RSC  array.  Current 
systems  leveraging  this  so-called  pupil-remapping 
architecture  include  the  FIRST  (see  e.g.  Perrin  et  al. 
(2006)  and  Huby,  E.  et  al.  (2013))  and  DRAGON¬ 
FLY  (Jovanovic  et  al.  2012)  instruments,  both  of 
which  were  designed  with  the  mission  of  imaging 
faint  objects  in  space  such  as  exoplanets.  A  phase  re¬ 
covery  algorithm  for  the  case  of  a  redundant  front-end  was 
proposed  by  Lacour  et  al.  (2007).  The  algorithmic  frame¬ 
work  developed  in  this  paper  is  directly  applicable  to  this 
architecture. 

While  the  classical  RSC  algorithms  cited  above 
have  been  applied  successfully  in  both  the  optical 
and  radio  realms,  their  applicability  is  limited  in 
the  former  relative  to  the  latter.  Observation  at  ra¬ 
dio  wavelengths  is  typically  limited  by  background 
noise  rather  than  by  shot  noise,  which  typically  lim¬ 
its  performance  at  optical  wavelengths.  As  a  result, 


along  a  given  triangle  of  spatial  frequencies.  We  will  use 
the  term  bispectra  to  refer  to  particular  samples  of  this 
function. 
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the  relative  strength  of  the  image  Fourier  compo¬ 
nents  (i.e.  the  visibility  function)  declines  more  grad¬ 
ually  with  increasing  baseline  length  in  the  former 
relative  to  the  latter  (Buscher  2015).  This  in  turn 
generally  makes  fringe  tracking  a  viable  option  for 
stabilizing  the  phase  measured  by  each  baseline  in  a 
radio  telescope  array,  and  in  such  cases  RSC-based 
inference  techniques  can  operate  directly  on  the  sta¬ 
bilized  Fourier  phasors.  We  will  refer  to  such  tech¬ 
niques  as  direct  techniques.  Direct  techniques  have  been 
implemented  as  the  calibration  scheme  of  choice  for  several 
radio  interferometers:  the  Donald  C.  Backer  Precision  Array 
for  Probing  the  Epoch  of  Reionization  (PAPER)  in  South 
Africa  (see  Ali  et  al.  (2015)),  the  MIT  Epoch  of  Reion¬ 
ization  (MITeOR)  in  the  United  States  (see  Zheng  et  al. 
(2014)),  and  the  Ooty  Radio  Telescope  (ORT)  in  India  (see 
Marthi  &  Chengalur  (2014)).  The  rapid  decay  in  SNR 
with  baseline  length  commonly  observed  at  optical 
wavelengths,  on  the  other  hand,  entails  integration 
of  atmosphere-invariant  observables  (e.g.  the  bispec¬ 
tra)  over  multiple  atmospheric  coherence  times  in 
the  case  of  dim-target  observation.  While  the  ob¬ 
servables  may  differ,  there  are  important  parallels 
to  be  noted.  Firstly,  the  ability  of  RSC  techniques 
to  recover  the  Fourier  phase  with  no  prior  knowl¬ 
edge  about  the  target  is  retained  with  atmosphere- 
invariant  observables  as  we  will  see.  Secondly,  in 
direct  and  atmosphere-invariant  approaches,  there 
are  two  principal  algorithmic  approaches  for  phase 
recovery:  those  that  operate  on  the  phasor  of  the 
observed  complex  visibilities  (respectively,  the  bis¬ 
pectra),  and  those  that  operate  on  their  phase  (re¬ 
spectively,  the  closure  phase).  An  important  issue 
encountered  in  both  approaches  is  the  issue  of  in¬ 
teger  phase  ambiguities,  which  manifests  itself  as 
phase  wrapping  in  phase-based  approaches,  and  false 
minima  in  the  phasor-based  approaches.  Array  de¬ 
sign  considerations  based  on  the  notion  of  wrap- 
invariance  (Kurien  et  al.  2016),  as  well  as  algorithms 
from  lattice  theory  (Lannes  &:  Anterrieu  1999),  have 
been  used  to  understand  and  solve  these  ambigu¬ 
ity  issues  in  the  case  of  direct  observation.  Our 
case  for  the  reliability  of  our  methods  will  include 
an  extension  of  these  principles  to  inference  from 
at  mosphere-  invariants . 

In  this  paper,  we  develop  a  novel  algorithmic 
framework  for  imaging  dim,  complex  objects  using 
atmosphere-invariant  observables  formed  by  a  re¬ 
dundant  configuration  of  a  large  number  of  optical 
telescopes  or  sub-apertures.  As  such,  this  framework 
may  prove  useful  to  astronomers  processing  data 
from  existing  systems  like  those  mentioned  above, 
as  well  as  to  the  designers  of  new  large-scale  in¬ 
struments  such  as  the  future  Planet  Formation  Im¬ 
ager  (PFI)  (Kraus  et  al.  2014),  which  seeks  to  de¬ 
velop  an  understanding  of  the  processes  governing 
the  very  assembly  of  planetary  systems  from  inter¬ 
stellar  dust.  To  the  best  our  knowledge,  work  by 
Lannes  &  Anterrieu  (1999)  was  the  first  to  explore 
the  use  of  such  atmosphere-invariant  observables  in 
the  explicit  context  of  RSC-based  phase  recovery. 
In  this  work,  the  authors  develop  a  pseudo-inverse  estima¬ 


tor  which  recovers  the  Fourier  phases  from  the  measured 
baseline  phases  via  an  intermediate  computation  of  the  clo¬ 
sure  phases.  The  computation  of  the  necessary  covariance 
matrices  for  optimal  estimation  within  this  formulation  is 
omitted.  Since  it  is  in  fact  the  closure  phases  which  are  the 
direct  observables  in  practical  low-flux  scenarios,  we  present 
instead  an  approach  working  solely  with  the  closure  phases 
and  compute  the  necessary  covariances  in  straightforward 
manner.  Lannes  (2003)  later  proposed  an  alternate  estima¬ 
tor  based  on  computation  of  the  Smith  Normal  Form  of  the 
matrix  mapping  Fourier  phases  to  closure  phases.  This  esti¬ 
mator  has  the  advantage  of  completely  eliminating  the  effect 
of  phase  wrapping  in  the  closure  measurements  when  appro¬ 
priate  routines  from  lattice  theory  are  used  to  pre-process 
the  phase  measurements.  In  contrast,  our  estimator  reliably 
uses  standard,  fast  linear  estimation  techniques,  thereby  ob¬ 
viating  computation  of  the  Smith  Normal  Form  at  the  pos¬ 
sible  expense  of  an  extra  shift  in  the  recovered  image.  Such 
image  shifts  are  anticipated  to  be  of  negligible  importance 
in  most  cases  of  practical  interest. 

An  important  novel  aspect  of  our  framework  is 
its  generalization  of  the  bispectrum  observable  for 
the  purposes  of  increased  sensitivity.  Though  the  bis¬ 
pectrum  and  its  associated  closure  phase  have  been  standard 
interferometric  observables  for  decades  after  the  latter 
was  first  proposed  by  Jennison  (1958),  it  is  not  difficult 
to  imagine  situations  in  which  use  of  these  observables  un¬ 
necessarily  limits  reconstruction  performance.  In  Figure  1, 
we  illustrate  this  idea  with  an  example.  We  label  each  base¬ 
line  of  an  interferometric  array  with  a  visibility,  which  is  an 
indicator  of  the  strength  of  its  Fourier  component  relative  to 
the  overall  brightness  of  the  image.  In  accordance  with  the 
power-law  decay  of  intensity  with  spatial-frequency  modu¬ 
lus  in  an  overwhelmingly-large  fraction  of  natural  images 
(see,  e.g.  Ruderman  (1994)),  we  consider  a  visibility  distri¬ 
bution  which  drops  sharply  with  baseline  length.  In  stan¬ 
dard  bispectral  imaging,  the  high-spatial-frequency  phase 
information  associated  with  long  baseline  b*  is  recovered  via 
forming  a  closure  triangle,  for  example,  with  a  short  base¬ 
line  along  with  another  long  baseline  ( Traditional ).  As  we 
will  see,  for  a  certain  range  of  photon  fluxes,  the  SNR  of  the 
bispectrum  is  roughly  proportional  to  the  sum  of  the  recip¬ 
rocals  of  the  squared  visibilities  of  the  associated  baselines 
(Kulkarni  et  al.  1991),  i.e.: 


SNR  »  n 


3  \ 


-1 


(2) 


where  h  is  the  number  of  photoelectrons  (pe)  received  per 
interference  fringe  per  exposure  frame.  As  we  will  show  in 
this  paper,  this  SNR  model  extends  in  the  natural  way  to 
the  SNR  of  higher-order  observables  like  that  shown  on  the 
right  ( Generalized ).  In  contrast  with  the  Traditional  ob¬ 
servable,  Generalized  observable  utilizes  only  short  (high- 
visibility)  baselines  to  close  the  long  baseline,  resulting  in  a 
near-doubling  of  the  resulting  SNR.  We  will  henceforth  refer 
to  the  higher-order  generalization  of  the  bispectrum  triple 
product  as  the  n-spectrum,  and  its  phase  as  the  generalized 
closure  phase.  The  statistics  for  these  new  observables  will 
be  derived  in  the  course  of  this  paper. 

Generalization  of  the  closure  phase  is  not  a  novel  con- 
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Figure  1.  Generalizing  the  phase  closure  concept  .  SNRs  given 
assume  each  aperture  contributes  n  =  2e3  photoelectrons  to  each 
fringe  in  pairwise  combination. 


cept  in  the  astronomical  community.  Lannes  (1991)  de¬ 
fined  the  general  class  of  matrices  which  annihilate 
the  subspace  comprising  all  possible  OPD  vectors. 
Recently  Martinache  (2010)  proposed  one  conve¬ 
nient  construction  of  such  matrices  as  the  orthog¬ 
onal  basis  for  the  left-nullspace  of  the  matrix  A  in 
Equation  (1).  The  vectors  of  this  basis  are  assembled  as 
the  rows  of  a  matrix  K.  When  this  matrix  is  applied  to  the 
measurements  /3,  the  atmosphere-induced  (/>  contributions 
clearly  cancel,  leaving  only  linear  combinations  of  the  object 
phases  in  0.  These  linear  combinations,  known  as  kernel 
phases ,  can  then  be  used  to  reconstruct  the  Fourier  phase 
via  direct  application  of  the  pseudo-inverse  of  K  (Martinache 
2014),  or  via  regularized  inversion  using  prior  image  models 
(Ireland  2013).  While  kernel  phases  can  serve  well  as 
observables  for  single-exposure  mapping,  they,  like 
their  closure  phase  counterparts,  may  fail  to  aver¬ 
age  reliably.  Indeed,  except  for  the  case  in  which 
the  kernel  phase  is  confined  to  a  narrow  region  in 
the  phase  domain  away  from  ±n  (e.g.  in  a  scenario 
in  which  the  per- frame  flux  is  high  and/or  the  target 
is  highly  compact),  the  inherent  2^-wrapping  of  the 
phase  renders  the  the  averaged  kernel  phase  a  poor 
estimator  of  the  true  kernel  phase  in  general  \  In 
this  paper  we  consider  the  scenario  of  observing  a 
complex  target  with  a  long-baseline  interferometer, 
in  which  the  histogram  of  measured  baseline  phases 
(and  hence  their  linear  combinations)  will  approach 
a  uniform  distribution  at  a  suitably-low  visibility. 
In  such  cases,  it  is  necessary  to  integrate  the  cor¬ 
responding  atmosphere-invariant  phasor  quantities, 
i.e.  the  n-spectra  defined  above,  and  then  extract 
their  phase  as  a  reliable  observable. 

Assuming  the  Fourier  phases  can  be  estimated  accu¬ 
rately  from  the  atmosphere-invariant  observables,  image 
reconstruction  essentially  reduces  to  the  same  frequency- 
undersampled  problem  encountered  in  a  plethora  of  other 
fields  including  medical  imaging  as  well  as  radio  astronomy. 
The  fundamental  problem  in  the  latter  is  to  regularize  the 


3  The  reader  is  directed  to  Buscher  (2015)  for  a  more  detailed 
discussion  of  this  issue. 


ill-posed  reconstruction  of  Nrespi  resolution  elements 
from  m  <K  Nrespi  spatial  frequency  measurements  corre¬ 
sponding  to  the  available  baseline  pairs  in  the  telescope  ar¬ 
ray.  Traditionally  approaches  based  on  the  so-called  CLEAN 
algorithm  due  to  Hogbom  (1974),  which  implicitly  fit  the 
measurements  to  an  image  model  consisting  of  a  sparse  col¬ 
lection  of  point-like  sources,  have  been  used  to  solve  this 
problem.  Given  the  success  of  this  relatively-simple  algo¬ 
rithm,  it  is  hardly  surprising  that  recently-developed  tech¬ 
niques  leveraging  sophisticated  sparse  image  models  (see  e.g. 
Wiaux  et  al.  (2009),  Baron  et  al.  (2010),  Kurien  et  al. 
(2014))  have  shown  promise.  In  fact  such  algorithms  be¬ 
long  to  a  burgeoning  family  of  cross-disciplinary  techniques, 
which  are  based  on  seminal  work  within  the  last  decade 
(Donoho  2006),  (Candes  et  al.  2006),  in  the  field  of  com¬ 
pressed  sensing. 

In  this  paper,  we  develop  a  comprehensive  algorithmic 
framework  for  image  reconstruction  based  on  RSC  tech¬ 
niques  and  compressed  sensing,  and  analyze  its  performance 
for  the  two  popular  beam-combining  architectures  used  in 
optical  interferometry  ( pairwise ,  and  Fizeau).  The  paper  is 
organized  as  follows.  In  Section  II,  we  derive  an  SNR  model 
for  the  11-spectrum,  as  well  as  for  the  covariance  amongst  dis¬ 
tinct  n-spectra.  The  analysis  begins  with  the  pairwise  case 
for  simplicity,  and  mathematical  arguments  for  an  extension 
to  the  Fizeau  case  are  presented  at  the  end  of  the  section 
and  in  the  Appendix.  In  Section  III,  we  describe  a  system¬ 
atic,  integer-least-squares  approach  for  reconstructing  im¬ 
agery  from  object  visibilities  and  generalized  closure  phases 
(GC).  We  draw  on  the  notion  of  a  minimum  cycle  basis  to 
select  a  minimum  set  of  linearly-independent  GC’s  (i.e.  a 
basis)  of  minimum  variance  which  spans  the  subspace  of  all 
closures.  Leveraging  techniques  first  suggested  by  Lannes 
&  Anterrieu  (1999)  and  the  notion  of  wrap-invariant  mea¬ 
surement  mappings  introduced  in  Kurien  et  al.  (2016),  we 
apply  algorithms  from  lattice  theory  to  unwrap  the  general¬ 
ized  closures.  We  then  quantify  the  gain  in  theoretical  per¬ 
formance  afforded  by  generalizing  the  classical  third-order 
atmosphere-invariant  observables  to  higher-order  cycles,  and 
corroborate  this  theoretical  performance  with  simulation  re¬ 
sults.  In  this  simulation,  we  apply  our  algorithm  to  recon¬ 
struct  a  dim,  structured  object  from  generalized  closures. 
For  the  final  image  reconstruction,  we  employ  a  sparse  re¬ 
covery  algorithm  which  enforces  smoothness  in  the  image 
domain  by  minimizing  the  aforementioned  total  varia¬ 
tion  metric. 

2  PRELIMINARIES 

2.1  Interferometric  Architectures 

In  this  paper  we  consider  two  popular  beam  combination 
architectures  in  use  in  optical  interferometry:  the  pairwise 
combination  scheme,  and  the  Fizeau  combination  scheme. 
The  two  schemes  are  illustrated  in  Figure  2  below.  Suppose 
we  have  a  telescope  array  consisting  of  Nap  apertures.  In 
the  pairwise  scheme,  light  from  each  aperture  is  split  Nap 
ways  and  combined  with  each  other  aperture.  Hence  each 
of  the  focal  planes  receives  a  small  fraction  of  the 

total  signal  N2" — j-  photons,  and  thus  incurs  the  shot  noise 
from  only  two  apertures.  In  the  Fizeau  scheme,  a  single  focal 
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Figure  2.  The  two  popular  beam  combination  schemes  in  optical 
interferometry 
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plane  receives  all  of  the  light  collected  by  the  array  (i.e.  Napn 
photons),  and  incurs  the  shot  noise  from  all  Nap  apertures. 

In  this  paper  we  analyze  the  performance  of  generalized- 
closure-based  imaging  for  both  architectures  in  terms  of 
mean-squared  error  (MSE)  of  the  phase  estimates.  We  begin 
with  the  pairwise  architecture  since  it  is  significantly  eas¬ 
ier  to  analyze.  We  then  present  a  series  of  approximations 
which  accurately  describe  the  algorithm’s  performance  in 
the  Fizeau  case  while  keeping  the  tediousness  of  the  math¬ 
ematics  at  a  manageable  level.  In  both  cases,  the  analysis 
is  based  on  the  computation  of  the  moments  of  the  Poisson 
distribution;  in  the  pairwise  case,  the  first  two  moments  are 
required,  whereas  in  the  Fizeau  case,  the  first  t  moments 
are  needed,  where  t  is  the  order  of  the  generalized 
closure. 


2.2  Fringe  Noise  Model  for  Pairwise  Beam 
Combiner 

Given  Nap  apertures,  recall  that  the  Fourier  phase  and  am¬ 
plitude  of  the  object  is  encoded  in  a  series  of  interfer¬ 

ence  fringes  observed  on  a  focal  plane.  In  the  pairwise  case, 
each  fringe  is  measured  on  a  separate  set  of  detectors  so 
that  the  Poisson  shot  noise  incurred  is  independent  for  each 
fringe.  Suppose  each  aperture  receives  n  photons  and  that 
this  light  is  split  evenly  (e.g.  by  a  beam-splitter)  before  be¬ 
ing  combined  with  the  light  from  the  other  apertures.  Hence 
each  aperture  contributes  h  =  ,,  -  photons  to  each  fringe. 
The  expectation  of  number  of  photons  q  at  a  given  detector 
k  for  a  given  fringe  is  given  by: 


To  obtain  the  bispectrum  SNR,  it  will  be  useful  to  com¬ 
pute  the  power  in  the  fringe,  i.e.  expectation  of  the  quantity 
* 

zz  ■■ 


( zz *} 


iNv  —  l  Nrj-l 


E  2  q(k)q(k')e-^ke^k’ 

k=0  k'= 0 


This  is  equivalent  to: 


(5) 


IVp-1  iVp-1 

<zz*>  =  E  E  m)q{k'))e-^ke^k'  (6) 

k=0  k'= 0 

Using  the  first  two  moments  of  the  Poisson  distribution, 
we  can  write: 


(q(k)q(k'))  =  < q(k))(q(k '))  +  Skk,(q(k ))  (7) 

Substituting  into  Equation  (6),  we  obtain: 

Np- 1  Np- 1  Np- 1 

<zz*>  =  E  ^me-i(ok  E  (9(k'))ei(ok'  +  £  <$(*)>  (8) 

k= 0  k'= 0  k= 0 

Simplifying  we  obtain: 

(zz*)  =  y2fr  +  2  h  (9) 


2.3  N-Spectrum  Covariance  and  SNR  Models  for 
the  Pairwise  Architecture 

Let  us  define  a  generalized  N-spectrum  G  as  the  product  of 
complex  fringe  phasors  along  an  f-edge  cycle.  Given  vari¬ 
ances  °~'Pe(G)  and  cr|  (G)  of  real  and  imaginary  compo¬ 
nents,  respectively,  we  define  the  pseudo-variance  of  the 
product  as: 


VpairwiseiG)  :=  cr2Re  +  =  < GG *)  -  (G)(G*) 

or: 

Vnni 


.(G)  = 


[>•<>  —  El  <^i><z*> 


(10) 


(11) 


r=i 


i=i 


V, 


pairwise 


(G)  = 


f]  («27?  +  2 h)  -  h2t  P|  yr 


(12) 


i=l 


i=l 


In  the  specific  case  of  the  classic  triple  product  (i.e.  the 
bispectrum ),  this  expands  to  (Kulkarni  et  al.  1991): 


(q(k))  =  — (1  +  ycos((x)k  +  6  +  (f>))  (3) 

where  Np  is  the  number  of  pixels  over  which  these  photons 
are  spread,  y  is  the  visibility  of  the  fringe  (which  takes  it 
value  between  0  and  1),  a>  is  the  fringe  frequency,  6  is  the 
Fourier  phase,  and  < j>  is  the  atmosphere  phase. 

To  retrieve  the  object  phase  and  amplitude  from  its 
fringe  encoding,  we  simply  take  the  Fourier  transform  and 
evaluate  it  at  the  fringe  frequency: 

Np- 1 

z=  E  <l^i0ik  (4) 

k= 0 


VpairwiseiG )  =  2n5  (y2yl+yjyj+y2y2)+4n4(y2+yl+y2)+Sn3 

(13) 


The  SNR  Spairwise  ■=  ^==  is  then  given  by: 

^  3 

C  _  717273^ 2 

^ pairwise  ~  , - 

yj»2(y \y\  +  yjyj  +  737?)  +  2n(y2  +  y\  +  y\)  +  4 

(14) 

To  gain  insight  into  the  limiting  behavior  of  the  SNR, 
we  make  two  approximations:  A  low-SNR  approximation  can 
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be  found  by  considering  only  the  constant  term  in  the  de¬ 
nominator,  which  will  dominate  when  h  and/or  the  y,-  are 
small. 


^  pairwise, low 


3 

717273" 2 
2 


(15) 


A  high-SNR  approximation  can  be  found  by  consider¬ 
ing  the  highest-order  term  in  the  denominator,  which  will 
dominate  when  h  is  large  and/or  the  y,-  are  large. 


3 

o  7l7273«2 

*^pairwise,high  ~  i - 

^2(r?  y\  +  y\y\  +  y\y\) 


(16) 


By  averaging  the  fringe  phasors  over  a  sufficient  number 
of  frames  Nf,  we  can  build  the  SNR  to  a  level  at  which  the 
bispectrum  phase  (i.e.  the  closure )  phase  can  be  measured. 
For  phasor  SNRs  sufficiently  greater  than  1,  the  following 
approximation  holds  for  the  variance  on  the  closure  phase 
Q el¬ 


and  G2  =  Re[G2\  +  jIm[G2 ],  and  £[G2l  =  g2-  Let  1  and  J  de¬ 
note  the  set  of  baselines  in  n-spectra  Gi  and  Gi,  respectively. 
Then  we  have: 

rrfG,.  G*)  =  (U  z,  f]  ^  ~  ft  <^>  II  <*?  ^ 

'iel  jej  1  iel  jej 

We  can  then  factor  out  the  factors  common  to  Gi  and 
Gi  from  those  which  are  distinct  to  obtain: 


o-fGpG*) 


n  n  n  ^  n 

.. kelnj  kelnj  lel\j  mej\l 

(25) 


Making  the  appropriate  substitutions  from  Section  1.2, 
we  obtain: 


del 


NfS2 


(17) 


We  then  obtain  the  following  approximations  for  the 
closure  phase  variances.  At  low  SNR,  we  have: 


Gci,low 


Taking  the  logarithm  we  find  that  the  variance  decou¬ 


ples: 

log  ~  ~2  log  71-2  log  72-2  log  73-3  log  h  +  C  (19) 

where  C  is  a  constant  independent  of  the  flux  level  h  and 
the  visibilities  7,-. 

At  high  SNR,  we  have: 


ecUhigh 


1  ,  1  1  1  . 
nNf  y\  y\  75 


(20) 


Generalizing  the  results  above  to  n-spectra,  we  obtain 
for  low  SNR, 


tr(G1,G*)  =  e2^'(ZG>-ZG2> 


FI  («27fc  +  2n6Zlk,Z2k) 
.kelnj 


n h 

kelnj 


x  n fiyt  n  "7m 

lel\j  mej\l 


where  SZlk?Z2k  =  1  if  the  common  measurements  are  the 
same,  and  zero  if  they  are  conjugates  of  each  other.  Similarly, 


o-(Gi.G2)  =  e^A^Gi+ZGi) 


|  I  Gi2yl+2h5Zl  k,Z2k 

.kelnj 


n r% 

kelnj 


x  n  "r/  n  "rm  *-27) 

lel\J  meJ\I 


where  6Zi  ktZ2  k  =  1  if  the  common  measurements  are  conju¬ 
gates  of  each  other,  and  zero  if  they  are  the  same. 

We  can  now  compute  real  and  imaginary  n-spectra  co- 
variance  as: 


0cl,loW  ~  Nfht  nj=1  7? 

and  for  high  SNR, 


(21) 


2  ~  1  1 
aGci,high  ~  Jjyj.  Z J  ^2 


(22) 


where  t  is  the  order  of  the  generalized  closure.  Note  that  if 
we  take  the  logarithm  of  Equation  (23),  we  obtain: 


cr(Re[Gil  Re[G2])  =  \_Re[Cov{Gh  G*)  +  Cov(Gh  G2)]  (28) 

o-(Im[GilRe[G2])  =  \lm[Cov(GhG*2)  +  Cov(GhG2)]  (29) 

o-(Re[GilIm[G2])  =  -Im[-Cov(GhG~)  +  Cov(GhG2)]  (30) 
2  z 


lQg  alcl,loW  «  f  log  -  -  2  ^  log  7i  -  log  Nf  (23) 

1  =  1 

2.4  Covariance  Matrix  of  the  Generalized  Closure 
Phases  for  the  Pairwise  Architecture 

In  this  section  we  generalize  results  of  Kulkarni  et  al.  (1991) 
to  obtain  expressions  for  the  n-spectra  covariance  matri¬ 
ces.  We  will  first  compute  the  covariance  between  two  n- 
spectrum  phasors  Gi  =  Re[G\]  +  jIm[G\\,  with  £T[G  1  ]  =  gi, 


<r(Re[Gi\,  Im[G2])  =  -Re[-Cov(Gh  G~)  -  Cov(Gh  G2)]  (31) 
2  z 

Given  the  n-spectra  covariance  expression  above,  we  can 
now  approximate  the  covariance  between  two  generalized 
closure  phases  using  Taylor  series  approximation  of  the  mo¬ 
ments,  which  is  a  technique  known  as  the  delta  method  (see, 
e.g.  Casella  &  Berger  (2002)).  First  we  define  the  angle  func¬ 
tion  of  a  phasor  as  6(x,  y)  =  arctan  L .  The  corresponding  gen¬ 
eralized  closure  phases  are  given  by:  =  6{Im[G\\,  Re[G\]) 

and  ©2  =  0(Im\_G2\,Re\G2\).  Applying  the  delta  method 
yields  the  following  approximation: 
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-<®..®2)-(^|g,)(5E)-(*«[G,l,MG2l)  (32) 

+(f!U(EJ'('”[G>ls'IC2l)  (33> 
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3  FOURIER  PHASE  RECOVERY  USING  THE 
N-SPECTRA 


Evaluation  of  the  derivatives  in  the  expression  above 
yields: 


cr(©l,©2)  ~ 


j  Re[go\  \ 

Re[g  ll2  +  Im[gi\2) 

\  A<?[g2]2  +  Im[g2]2) 
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|  /  -Im[g2] 
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1 7?e[g2J2  +  Im[g2]2 
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l  ~Im[g2\ 
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<x(I?e[Gi],  7m[G2]) 
Jo-(/m[Gi],  Re[G2]) 
Jo-(/m[Gi],  7m[G2]) 
)  cr(/fe[Gi],  Re[G2}) 


(36) 


2.5  Covariance  Matrix  Approximations  for  the 
Fizeau  Architecture 


In  this  section,  we  briefly  review  the  integer  least-squares  ap¬ 
proach  for  RSC  phase  recovery  (see  e.g.  Lannes  &  Anterrieu 
(1999),  Kurien  et  al.  (2016)  for  more  comprehensive  treat¬ 
ments),  and  then  describe  an  algorithm  for  selection  of  a  set 
of  generalized  closure  relations  of  minimum  total  variance. 


3.1  RSC  Phase  Recovery  with  Wrap-Invariant 
Patterns 

The  traditional  approach  to  RSC  reconstruction  operates  on 
the  measured  baseline  phases  (see  e.g.  Arnot  et  al.  (1985), 
Greenaway  (1994)).  To  illustrate  the  approach,  let  us  con¬ 
sider  an  interferometer  which  operates  at  a  wavelength  A 
with  two  apertures  (say,  i  and  j)  separated  by  a  vector  base¬ 
line  distance  of  by.  In  the  absence  of  any  optical  path  differ¬ 
ence,  the  interference  pattern  formed  by  these  two  apertures 
encodes  a  sample  of  the  object’s  Fourier  Transform  at  spatial 
frequency  -j- .  Let  the  true  Fourier  phase  (which  we  will  re¬ 
fer  to  as  object  phase),  measured  by  this  interference  pattern 
be  denoted  as  0y .  The  measured  phase  is  given  by: 


In  the  Fizeau  case,  Equation  (9)  becomes: 


Pi]  =  %  +  0/  _  <Pi  +  2ne 


(42) 


(zz*)  =  y2n2  +  Napn  (37) 

As  shown  in  Appendix  A,  we  can  approximate  the  vari¬ 
ance  of  a  n-spectrum  of  order  o  as: 


V Fizeau  (G) 


f]  (n2y2  +  Napn)  ~n2tfjy- 


1=1  1=1 
Hence  the  SNR  can  be  approximated  as: 


(38) 


•Sri 


V2«r  n'=1  r  i 


V[n -=i  +  Napn)]  -  «2r  n[=1  yj 


(39) 


Applying  the  same  approximations  in  Appendix  A  to 
the  covariance  yields  the  following  analogues  to  Equations 
(26)-(27) 


o-(GhG*2)  =  ^iCGrZG)) 


n  (»2yf  +  Napn6ZlkiZ2k)-  ]~~[  n2y2  ]"“[  ny,  ]““[  nym 
kelnj  kelnj  lel\j  mej\l 


(40) 


and, 


where  <t>j-<pi  is  the  optical  path  difference  between  apertures 
j  and  i,  and  e  is  unknown  phase  wrap  integer  arising  from 
the  fact  that  interferometric  phase  measurements  are  only 
known  modulo  In. 

Consider  an  interferometric  array  which  simultaneously 
makes  many  such  measurements  amongst  its  Nap  apertures. 
Suppose  that  of  the  array’s  (iV2p)  baselines,  d  of  them  are 
distinct.  Further  suppose  we  have  a  solution  set  (</>,-)  and 
{ Ofj }  for  these  equations.  Let  r,-  denote  the  vector  position 
of  the  i-th  aperture.  As  noted  by  several  authors  (see,  e.g. 
Wieringa  (1992)),  we  can  obtain  another  valid  solution  set 
simply  by  replacing  each  <pi  with  ^7  =  ^,  +0o  +  z'rii  and  each 
Ojj  with  0?j  =  djj  -  z  •  (r j  -  r ;),  for  arbitrary  </>q  and  z.  Since 
the  free  vector  z  is  a  two-parameter  vector  representing  the 
inherently-ambiguous  position  of  the  image  within  the  field- 
of-view  and  the  free  parameter  <f>o  is  simply  a  scalar  piston 
offset,  the  kernel  of  the  RSC  system  is  three-dimensional. 
Therefore  the  RSC  system  contains  d  unknown  distinct  ob¬ 
ject  phases  and  Nap  unknown  aperture  pistons,  and  is  rank- 
deficient  by  at  least  3.  This  implies  that  there  are  at  most 
(d  +  Nap  -3)  linear ly-independent  equations  in  the  RSC  sys¬ 
tem,  and  hence  at  most  Nap  -  3  redundant  relations  can  be 
linearly-independent.  We  will  assume  for  the  remainder  of 
the  paper  that  our  array  contains  Nap  -  3  independent  rela¬ 
tions.  Under  this  assumption,  we  can  then  solve  for  a 
particular  solution  of  this  system  by  arbitrarily  set¬ 
ting  the  value  of  two  object  phases  (whose  spatial 
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frequencies  are  not  co-linear)  and  one  piston  phase. 

This  particular  solution  will  then  differ  from  the  true  solu¬ 
tion  by  a  phase  ramp  in  the  Fourier  domain,  corresponding 
to  an  image  shift  in  the  spatial  domain. 

Using  M  to  denote  the  coefficients  of  the  equations,  we 
can  write  the  RSC  system  in  compact  form  as: 


M 


0 

$ 


/3  +  2ne 


(43) 


We  now  show  how  generalized  closures  eliminate  the 
0  contribution  to  the  model  in  Equation  (43)  and  thereby 
provide  a  means  for  Fourier  phase  recovery  in  sce¬ 
narios  where  fringe  tracking  capabilities  are  either 
unavailable  or  degraded  due  to  low  SNR.  First  we 
define  the  interferometric  graph  of  an  array  as  the  directed 
graph  whose  vertices  are  the  Nap  apertures  in  the  array  and 
edges  are  the  | p)  baselines  connecting  all  vertex  pairs.  We 
also  define  the  spectral  phase  space  K  as  the  span  of 
the  d  columns  associated  with  the  Fourier  phases,  or  equiv¬ 
alently  the  range  of  the  sub-matrix  Mg  formed  by  these 
columns.  The  piston,  or  aberration ,  phase  space  L  is  the 
span  of  the  Nap  columns  associated  with  the  piston  phases, 
or  equivalently  the  range  of  the  sub-matrix  formed  by 
these  columns.  □ 

Note  that  the  d  columns  spanning  K  are  linear ly- 
independent  by  virtue  of  having  non-zero  entries  in  mutually 
disjoint  sets  of  baseline  indices.  Hence  dim(K )  =  d.  The  sub¬ 
space  L  has  dimension  Nap  -  1  (Lannes  &  Anterrieu  1999); 
the  constant  vector  forms  the  one-dimensional  nullspace  of 
M*. 

We  generalize  the  notion  of  phase  closure  in  a 
structured  manner  by  defining  the  important  con¬ 
cept  of  the  cycle  space  of  the  interferometric  graph. 

Given  a  directed  graph  G  with  a  set  of  vertices  V  and  edges 
E,  the  cycle  space  of  the  graph  is  the  vector  space  of 
spanned  by  the  incidence  vectors  of  cycles  with  a  (clockwise, 
or  counter-clockwise)  orientation.  It  can  be  shown  that  the 
dimension  of  the  cycle  space  is  M  +  N  -  1,  where  M  =  \E\  and 
N  =  \V\. 

Suppose  we  have  a  basis  for  the  cycle-space  of  the  in¬ 
terferometric  graph.  Let  us  stack  the  elements  of  this  ba¬ 
sis  as  row  vectors  in  a  )  x  matrix  Cmc,  which 

maps  baseline  phase  measurements  to  generalized  closure 
phases.  Recall  from  the  discussion  of  n-spectra  in  Section 
1  that  closure  relations  eliminate  piston  differences  in  the 
measurements  so  that  Cmc  annihilates  the  subspace  L  (see 
e.g.  Lannes  (2003)),  i.e.  the  space  spanned  by  the  columns 
of  M  corresponding  to  (j>. 

In  Figure  3,  we  show  a  simple  four-aperture  interfero¬ 
metric  graph  and  one  possible  cycle  basis  for  the  graph.  For 
this  particular  example,  we  can  write: 


C 


me 


11  1  0  0  0 
0  0-1110 
10  0  011 


(44) 


where  the  column  (i.e.  baseline)  indexing  in  the  matrix  fol¬ 
lows  the  labeling  in  the  Figure.  Note  that  the  final  triangle 
cycle  (i.e.  the  one  containing  baselines  2,  4,  and  6,  and  rep¬ 
resented  by  the  row  vector  w*  :=  (0, 1, 0, 1, 0, -1)),  can  be 
represented  as  a  linear  combination  of  the  cycles  in  the  ba¬ 
sis:  w*  =  C[  +  C2  -  C3,  where  (c,- j  are  the  rows  of  Cmc 


Figure  3.  One  possible  cycle-basis  for  a  simple  interferometric 
graph.  Note  the  fourth  triangle  (i.e.  baseline  set  { 2,  4,  6 } )  can  be 
expressed  as  a  linear  combination  of  the  cycles  shown. 


While  for  purposes  of  simplicity  we  have  chosen  a  cycle 
basis  consisting  exclusively  of  three-baseline  cycles,  in  gen¬ 
eral  the  elements  of  a  cycle  basis  can  contain  any  number  of 
edges. 

As  the  closure  relations  defined  by  Cmc  annihi¬ 
late  the  subspace  L,  we  have: 


0 

=  [  Coc  0  1 

9  ' 

.  $  . 

$  . 

L  J 

(45) 


where  C „c  :=  CmcMg,  and  is  the  mapping  between  object 
phases  and  closure  phases. 

By  standard  principles  of  linear  algebra,  we 
know  that  the  ability  to  uniquely  recover  the  set 
of  object  phases  given  closure  phase  measurements 
amounts  to  the  injectivity  of  the  matrix  Coc,  which 
can  be  verified  by  examining  its  nullspace.  Let  Ar  v  and 
Ary  denote  the  vectors  containing  the  x-  and  v-coordinates 
of  the  baselines  in  an  array,  respectively.  If  the  array  is  a 
valid  RSC  array,  these  columns  form  a  basis  for  the  two- 
dimensional  nullspace  of  Coc.  Given  that  the  nullspace 
is  composed  of  linear  combinations  of  the  array  base¬ 
line  vectors,  which  are  in  turn  scaled  versions  of  the 
array’s  spatial  frequencies,  the  mapping  Coc  is  injec¬ 
tive  up  to  a  linear  frequency  ramp  (or,  equivalently, 
an  image  shift).  The  proof  is  given  in  Appendix  B. 

In  addition  to  ensuring  injectivity  of  the  closure 
mapping,  we  must  also  address  the  difficulty  aris¬ 
ing  from  the  fact  that  generalized  closure  phases  are 
only  known  In.  Let  us  define  yc./  as  the  observation  vec¬ 
tor  of  wrapped  generalized  closure  phases,  and  2 neci  the 
wrapping  vector  (with  integer  ec/).  Our  phase  measurement 
model  is  then  given  by: 


Coc9  =  ye  +  n  (46) 

where  ye  =  yci+2neci ,  and  n  is  the  measurement  noise,  which 
we  shall  assume  in  this  paper  is  predominantly  due  to  shot 
noise  on  the  focal  plane. 

We  can  now  formulate  the  well-posed  recovery 
of  the  Fourier  phases  as  the  following  generalized,  in¬ 
teger  least-squares  problem,  which  finds  a  vector  in 
the  range  of  the  closure  matrix  at  a  minimum  Ma- 
halanobis  distance  from  the  actual,  wrapped  mea¬ 
surement  vector: 

(>RSC  =  ar8minechg{ ye  ~  Coce)T£_1  (ye  -  C ocS)  (47) 

where  £  is  an  estimator  of  the  covariance  matrix  of  the  phase 
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measurements.  This  estimator  can  be  obtained  by  substitut¬ 
ing  estimates  for  the  object  visibilities  jy,  j  into  the  covari¬ 
ance  expressions  in  the  previous  section.  It  is  acknowledged 
that  this  least-squares  solution  is  neither  the  optimal  nor 
maximum-likelihood  solution  to  the  phase  inference  prob¬ 
lem.  For  one,  a  Gaussian  distribution  on  the  closure  phases 
is  implicitly  assumed  to  hold,  and  can  only  approximately 
hold  for  high-SNR  closure  phases.  It  is  therefore  assumed 
that  n-spectra  will  be  integrated  for  a  sufficiently  large  num¬ 
ber  of  frames  for  this  approximation  to  be  reliable  for  most, 
if  not  all,  of  the  closure  phases. 

Assuming  the  covariance  matrix  estimator  is  not  de¬ 
generate,  it  will  admit  a  Cholesky-decomposition  T  =  BB7} 
Equation  (47)  is  then  equivalent  to  searching  for  vectors  e 
which  minimize  the  projection  of  a  whitened  measurement 
B_1ye  onto  the  subspace  defined  by  ker(( B-1Coc)r).  Specif¬ 
ically  we  seek  to  minimize: 

/(e)  =  l|PsB_1(yc/  +  2neci)\\2  (48) 

where  Pi  is  a  matrix  representing  the  orthogonal  projection 
from  R"  onto  ker((B~1Coc)T). 

Letting  e'.;  =  -ec/,  we  can  rewrite  the  above  objective 
function  as: 

f(e')  =  ||PxB-1(ycz  -  27re'./)||2  =  UPfcB"^,  -  27rPsB- 1  1 12 

(49) 

This  optimization  problem  is  equivalent  to  the  so-called 
closest  vector  problem  (CVP)  in  the  theory  of  lattices.  The 
connection  of  interferometric  phase-unwrapping  with  the 
CVP  problem  was,  to  the  best  of  our  knowledge,  first  ex¬ 
plored  in  Lannes  &  Anterrieu  (1999).  In  this  paper,  the  au¬ 
thors  formulate  an  analogous  minimization  problem  for  the 
raw  baseline  measurements. 

We  will  define  a  lattice  L(Zn)  as  the  set  of  points  gen¬ 
erated  by  integer  combinations  of  the  column  vectors  of  a 
matrix  L.  Letting  P  =  P^B-1  and  factoring  2 n,  our  optimiza¬ 
tion  problem  then  amounts  to  the  following:  Find  the  lattice 
point  in  P(Z")  which  is  closest  to  Pyc//(2;r).  A  compact  rep¬ 
resentation  of  the  lattice  T  is  given  by: 

m<n—(d+N— 3)  i 

X  a/v<  I  Vfli  6  z  l  (50) 

where  (v,- )  are  linearly-independent  and  together  form  a  basis 
of  the  lattice.  Given  the  lattice  basis,  several  algorithms  ex¬ 
ist  for  finding  the  closest  lattice  point  to  a  specified  vector.  A 
popular  class  of  algorithms,  known  as  the  Sphere-Decoding 
algorithms,  are  efficient  searches  for  the  closest  lattice  point 
within  a  hypersphere  of  a  certain  radius  centered  on  the  in¬ 
put  vector  (see  e.g.  Agrell  et  al.  (2002)).  For  the  simulations 
in  this  paper,  we  instead  use  the  lower-complexity  Babai 
Nearest  Plane  (Babai-NP)  algorithm  (Babai  1986).  For  lat¬ 
tice  bases  which  are  nearly  orthogonal  (such  as  those  we  use 
for  our  simulations),  this  algorithm  offers  reliable,  albeit  not 
guaranteed,  performance  in  practice. 

Suppose  we  have  found  a  basis  for  the  lattice  P(Z"), 
and  we  have  solved  the  Closest  Vector  Problem  for  a  given 
closure  vector  yc/ .  Let  b*  be  the  output  of  the  Babai  Nearest 
Plane  Algorithm  -  i.e.  it  is  the  lattice  point  which  is  the 


closest  to  yc/.4  We  now  seek  to  solve  for  the  wrap  vector 
corresponding  to  this  lattice  point,  i.e.  we  seek  a  solution 
to: 

b*  =  Pe  (51) 

Note  that  P  is  a  projection  matrix  and  thus  not  full- 
rank,  and  therefore  there  will  be  infinitely-many  solutions 
to  this  equation.  Indeed  our  solution  y*.  is  correct  only  to 
within  an  integer  vector  eresjd  in  the  range  of  Coc.  Hence 
there  is  a  fundamental  ambiguity  that  needs  to  be  addressed. 
Once  the  measurement  vector  is  unwrapped,  there  are  var¬ 
ious  ways  to  solve  the  RSC  system  which  respect  the  pres¬ 
ence  of  a  residual  ambiguity  described  above.  Follow¬ 

ing  standard  least-squares  principles,  we  first  compute  via 
projection  the  vector  y*;  im(c)  in  /m( Coc)  closest  (in  Ma- 
halanobis  distance)  to  the  unwrapped  measurement  vector, 
i.e. 

yW)  =  B-1«-p2)By;z  (52) 

Per  the  first  method  (Lannes  2003),  we  can  compute 
the  Smith  Normal  Form  (SNF)  (Smith  1861)  of  Coc.  This 
matrix  decomposition  described  in  the  following  Theorem: 
Theorem  3.1  (Smith  Normal  Form)  (Smith  1861):  Let  A 
be  a  nonzero  m  x  n  integer  matrix  with  rank  r.  There  ex¬ 
ist  integer  and  unimodular  5  (and  thus  invertible)  matrices 
m  x  m  and  n  x  n  matrices  U  and  V  respectively  such  that 
the  matrix  product  D  =  UAV  is  a  diagonal  matrix  whose 
diagonal  entries  D,-,-  (the  so-called  elementary  divisors)  are 
non-zero  integers  for  i  <  r,  and  zero  for  i  >  r.  Moreover,  the 
matrices  U  {and  respectively,  the  matrix  V}  represent  the 
row  {and  column}  operations  which  zeroize  A  below  {and 
above}  the  diagonal.  □ 

Let  us  compute  the  Smith  Normal  Form  (SNF)  (U,  D,  V) 
of  our  matrix  Coc,  and  set  Uc  =  U-1,  =  D,  and  Yc  = 

v-1. 

Coc  =  UCDCVC  (53) 

A  valid  RSC  solution  can  then  be  obtained  by  evaluat¬ 
ing: 

dRSC=\-cl  D£U-ycUm(c)  (54) 

where  Djt,  denotes  the  pseudo-inverse  of  D. 

The  effect  of  the  residual  ambiguity  vector  tresid  on 
this  solution  requires  careful  consideration.  Note  that  the 
resulting  wrap-induced  error  in  this  solution  is  given  by: 

2?r eRSC  =  2xVc  Dcuc  eresid  (55) 

Since  the  matrices  and  are  unimodular,  they  are 

4  It  is  well-known  that  the  performance  of  the  Babai  algorithm 
and  other  CVP  algorithms  is  improved  when  the  lattice  basis  is 
reduced  to  one  which  has  short  and  nearly-orthogonal  vectors. 
We  used  the  LLL  algorithm  implementation  due  to  Zhou  (2014) 
in  our  simulations. 

5  A  unimodular  matrix  is  a  square  integer  matrix  whose  deter¬ 
minant  is  ±1 
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invertible  over  the  integers.  If  all  the  elementary  divisors  in 
Dc  are  equal  to  1,  then  will  also  be  integral,  and  hence 
eRSC  will  be  integral.  This  in  turn  guarantees  that  the  fi¬ 
nal  error  in  the  Fourier  phase  will  be  0  mod  2n.  i.e.  that 
the  RSC  solution  9rsc  is  wrap-invariant.  We  hence  refer  to 
the  mapping  Coc  as  a  wrap-invariant  closure  mapping.  In 
Kurien  et  al.  (2016),  it  was  observed  that  this  condition  can 
be  violated  by  RSC  patterns  which  are  not  wrap-invariant. 
Conversely,  we  have  observed  that  with  wrap-invariant  pat¬ 
terns,  the  closure  mappings  associated  with  the  cycle  bases 
that  we  will  consider  in  this  paper  (i.e.  fundamental  cycle 
bases  and  minimum  cycle  bases)  are  also  typically  wrap- 
invariant  in  practice. 

As  it  turns  out,  for  such  wrap-invariant  mappings,  the 
inverse  problem  in  Equation  (46)  can  be  solved  reliably  even 
without  the  Smith  Normal  Form;  it  is  sufficient  to  find  any 
unimodular  r-by-r  sub-matrix  C  of  Coc  and  solve  the  as¬ 
sociated  smaller  system  in  Equation  (56)  to  obtain  a  valid 
solution.  Recall  that  for  valid  RSC  arrays  Coc  will  have  r  +  2 
columns  and  hence  in  this  approach,  two  (non-collinear)  ob¬ 
ject  phases  are  implicitly  set  to  zero.  This  implicit  selection 
then  fixes  the  fundamentally-ambiguous  translation  of  the 
scene  discussed  at  the  beginning  of  the  section. 


Since  C  is  unimodular,  C-1  will  have  solely  integral 
entries  so  that  the  resulting  wrap  error  C_lerei!(/  =  0 
mod  2n.  Assuming  the  measured  phase  has  been  correctly 
unwrapped,  the  solution  in  Equation  (56)  amounts  to  a  so- 
called  basic  solution  of  our  generalized  least-squares  prob¬ 
lem.  We  then  have  the  following  expression  for  the  error 
covariance  matrix  for  the  estimator  9rsc  (Kay  1993). 

Zbasic  =  (Cr£-1C)-1  (57) 

The  basic  solutions  belong  to  the  countably-infinite  set 
of  solutions  to  Equation  (46).  It  turns  out,  somewhat  sur¬ 
prisingly,  that  any  particular  solution  within  this  set  can  be 
reliably  limited  to  a  mere  image  shift  for  patterns  that  are 
wrap-invariant.  The  complete  set  is  given  by: 

9cr  =  C+cycUm(c)  +  9q  (58) 

where  9q  is  any  vector  in  the  nullspace  of  Coc,  and  Cf}c  de¬ 
notes  the  pseudo-inverse  of  Coc.  The  pseudo-inverse  can  be 
computed  using  the  Singular- Value-Decomposition  (SVD)  of 
Coc,  which  is  given  by: 

Coc  =  (59) 

where  Uo-  and  Vo-  are  O'*  (^2  *)  and  orthogonal 

matrices,  respectively.  So-  is  a  ( 'Jxrf  diagonal  matrix 
with  r  non-zero  diagonal  entries  (the  so-called  singular  val¬ 
ues  of  Coc  ),  where  r  =  rank(Coc )  =  d  —  2. 

The  Moore-Penrose  pseudo-inverse  is  then  given  by 
(Bretscher  2001): 

Coc  =  U  (60) 

where  Ej.  is  a  diagonal  matrix  whose  r  non-zero  diagonal 


entries  are  the  reciprocals  of  the  corresponding  non-zero  en¬ 
tries  in  £o-- 

In  our  companion  paper  (Kurien  et  al.  2016),  we  show 
that  the  effect  of  wrapping  on  the  RSC  pseudo-inverse  so¬ 
lution  to  Equation  (43)  can  be  reliably  limited  to  an  image 
shift  for  wrap-invariant  patterns.  An  analogous  result  ap¬ 
plies  to  closure-based  RSC  imaging: 

Proposition  3.2:  If  a  pattern  is  wrap-invariant,  the  er¬ 
ror  induced  by  wrapping  of  the  (generalized)  closure  phases 
can  be  limited  to  an  image  shift. 

Proof:  See  Appendix  C  □ 

The  covariance  of  elements  in  the  pseudo-inverse  solu¬ 
tion  can  be  expressed  as  (Montgomery  et  al.  2006): 

Zpinv  =  (61) 

where  Vo-  is  obtained  by  omitting  the  final  2  columns  of  Vq- 
which  form  a  basis  for  the  nullspace  of  Coc. 

3.2  Selection  of  the  N-Spectra 

In  this  section,  we  describe  a  strategy  for  obtaining  a 
near-optimal  linearly-independent  set  of  generalized  closure 
phases.  This  strategy  is  founded  on  the  notion  of  minimum 
cycle  basis  from  graph  theory,  which  is  defined  as  fol¬ 
lows.  Let  each  edge  of  the  graph  be  assigned  a  positive 
weight  C£,  and  the  weight  of  a  cycle  be  defined  as  the  sum 
of  the  weights  of  its  constituent  cycles.  We  can  then  ask  for 
the  set  of  linearly-independent  oriented  cycles  which  spans 
the  cycle  space  and  has  minimum  total  weight.  We  call  such 
a  set  the  minimum  cycle  basis  of  a  directed  graph  G. 

The  canonical  method  for  finding  a  minimum 
cycle  basis  is  known  as  the  greedy  algorithm,  which 
is  defined  as  follows.  Consider  a  vector  space  V,  a  set  of 
vectors  which  span  V,  and  a  corresponding  set  of  weights  for 
each  of  these  vectors.  A  basis  of  minimum  total  weight  can 
be  found  by  first  sorting  the  vectors  in  decreasing  order  of 
weight,  and  then  selecting  vectors  for  the  basis  in  this  or¬ 
der  if  and  only  if  they  are  linearly-independent  of  previous 
selections.  It  can  be  shown  that  this  strategy  is  guar¬ 
anteed  to  output  a  minimum  cycle  basis.  (Gross  & 
Yellen  2006) 

We  can  find  a  set  of  linearly-independent  closure  rela¬ 
tions  of  minimum  total  variance  by  employing  the  greedy 
algorithm  on  the  set  S  of  all  possible  closure  relations.  How¬ 
ever,  this  set  is  equivalent  to  the  set  of  all  possible  cyclic 
permutations,  which  becomes  enormous  for  arrays  of  large 
size.  Fortunately,  in  an  extension  of  a  result  due  to  Hor¬ 
ton  (1987)  for  undirected  graphs,  Liebchen  &  Rizzi  (2005) 
showed  that  the  elements  of  the  minimum  cycle  basis  of  a 
directed  graph  could  be  found  among  a  special  (and  much 
smaller)  subset  of  5  known  as  the  Horton  cycles. 

To  define  the  Horton  cycle,  we  must  first  in¬ 
troduce  the  concept  of  shortest  path  tree  for  those 
readers  who  may  not  familiar  with  it.  The  shortest 
path  tree  of  a  graph  G  with  respect  to  a  node  A  is  the  span¬ 
ning  tree  which  connects  node  A  to  each  other  node  in  G  via 
a  path  of  minimum  weight.  Figure  4  gives  a  simple  example 
of  a  shortest  path  tree  for  an  interferometric  graph  in  which 
the  top-most  node  is  the  root.  The  construction  of  shortest 
path  trees  is  a  well-studied  problem  in  graph  theory  (see 
e.g.  Gross  &  Yellen  (2006))  which  is  typically  solved  using  a 
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Figure  4.  A  shortest  path  tree 


- Shortest-Path  Tree 


shortest-path  routine  such  as  Dijkstra’s  Algorithm  (Dijkstra 
1959). 

Given  the  shortest  path  tree  T  originating  from  an  ar¬ 
bitrary  node  A,  a  Horton  cycle  is  a  cycle  formed  by  connect¬ 
ing  the  endpoints  of  any  two  branches  of  7.  As  it  turns  out, 
all  elements  of  the  minimum-cycle-basis  are  Horton  cycles. 
Therefore  in  searching  for  the  minimum-cycle-basis,  it  suf¬ 
fices  to  search  the  Horton  cycles.  A  proof  of  this  result 
is  given  in  Liebchen  &  Rizzi  (2005). 

We  now  leverage  this  fundamental  result  from 
graph  theory  to  develop  a  practical  method  for  se¬ 
lecting  a  robust  set  of  atmosphere-invariants.  Specif¬ 
ically,  we  seek  a  decomposition  of  the  weight  of  a  cycle  into 
the  weights  of  its  constituent  edges  (i.e.  baselines).  As  the 
following  Propositions  show,  at  the  extreme  low-  and  high- 
SNR  limits,  the  n-spectrum  variance  does  indeed  allow  such 
a  decoupling. 

Proposition  3.3  ( Low-SNR  Decoupling  Approxi¬ 
mation):  In  the  low-SNR  limit,  the  logarithm  of  the  total 
variance  of  any  cycle  is  given  (up  to  a  global  additive  con¬ 
stant)  by  the  sum  of  the  costs  assigned  to  each  baseline  in 
the  cycle,  where  the  cost  of  a  baseline  k  with  visibility  jk  is 
given  by:  ck  =  log  §  -  2  log  yk . 

Proof:  This  can  be  seen  by  re-writing  the  low-SNR  ap¬ 
proximation  in  Equation  (23)  as: 

lo8  alchloW  x  Z  log  2  log  V*  +  C  (62) 

1  =  1 

where  C  =  log  Nf  is  a  constant  with  respect  to  both  visibility 
and  flux  and  therefore  irrelevant  for  the  purposes  of  closure 
selection.  □ 

Proposition  3.4  ( High-SNR  Decoupling  Approxi¬ 
mation):  In  the  high-SNR  limit,  the  total  variance  of  any 
cycle  is  given  (up  to  a  global  scaling  factor)  by  the  sum  of 
the  costs  assigned  to  each  baseline  in  the  cycle,  where  the 
cost  of  a  baseline  i  with  visibility  y;  is  given  by:  c,-  =  - j. 

'  i 

Proof:  This  follows  directly  from  the  high-SNR  approxima¬ 
tion  in  Equation  (22).  □ 

Based  on  the  Propositions  above,  we  now  have  a  way  to 
construct  a  set  of  linearly-independent  generalized  closures 
of  minimum  variance  at  the  low-  and  high-  SNR  extremes. 
As  demonstrated  in  the  next  section,  even  at  moderate  flux- 
levels  where  these  approximations  lack  precision,  this  new 
closure  selection  method  can  yield  significant  improvements 
over  traditional  closure  imaging  restricted  to  baseline  tri¬ 
angles.  To  estimate  the  scaling  of  the  computational  cost 
of  the  method,  recall  that  the  greedy  algorithm  is  run  on 
the  set  of  Horton  cycles.  Since  there  are  Nap  shortest-path 


trees,  each  rooted  at  a  distinct  aperture,  and  m  -  ( Nap  -  1) 
cycles  in  each  of  these  trees,  the  total  number  of  Horton 
cycles  is  0(mNap).  Since  m  =  0(N^p),  we  see  that  the  al¬ 
gorithm  is  0(Nap),  which  is  the  same  complexity  as  if  the 
greedy  algorithm  were  run  on  an  exhaustive  listing  of  all 
three-baseline  closures. 


3.3  Complete  Algorithm 

We  are  now  in  a  position  to  synthesize  the  techniques 
developed  in  previous  sections  into  an  algorithm  for  RSC 
imaging  using  generalized  closures.  Recall  that  closure  selec¬ 
tion  relies  on  estimates  of  the  baseline  visibilities  {y,  j.  These 
estimates  can  be  derived  for  the  pairwise  case  by  solving 
Equation  (9)  for  y  to  yield: 


7  i  = 


V(zz*>  -  2 n 

n 


(63) 


Similarly,  for  the  Fizeau  case,  we  can  solve  Equation 
(37)  to  obtain: 


7i  = 


fizz*)  -  Napn 

n 


(64) 


We  can  now  combine  corresponding  estimates 
of  the  Fourier  phases  and  visibilities  to  form  Fourier 
phasor  estimates  y,  which  are  known  as  complex  vis¬ 
ibilities.  To  estimate  the  image  coefficients  x  from  these 
complex  visibilities,  we  must  apply  a  deconvolution  algo¬ 
rithm.  Based  on  the  demonstrated  effectiveness  of  recent 
compressed  sensing  techniques,  we  choose  the  sparse  re¬ 
covery  technique  known  as  Total  Variation  Mini¬ 
mization  (Rudin  et  al.  1992)  to  perform  this  deconvolu¬ 
tion.  Namely  our  image  estimate  is  given  by: 


x  =  argmin^  Hallyy  subject  to:  ||Fcr  -  y||2  <  e  (65) 

where  F  is  a  partial  Fourier  matrix  whose  rows  are  vector¬ 
ized  representations  of  the  2D  sinusoids  associated  with  the 
array’s  measured  spatial  frequencies.  For  the  image  re¬ 
constructions  conducted  below,  we  used  the  NESTA 
software  package  (Becker  et  al.  2011)  to  perform  this 
regularization. 

Algorithm  1  provides  the  complete  set  of  steps  we  have 
discussed. 


4  ALGORITHM  PERFORMANCE 

In  this  Section,  we  present  the  results  of  application  of  Al¬ 
gorithm  1  to  a  simulated  scenario  of  imaging  a  structured 
object  space. 


4.1  Sensitivity  Limits 

A  useful  benchmark  for  our  results  is  provided  by  the 
Cramer-Rao  Lower  Bound  (CRLB)  for  interferometric  phase 
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Algorithm  1  RSC  Image  Reconstruction  Algorithm 
1.  Estimate  visibilities  y,-  (c.f.  Equations  (63)-(64)) 

2  Select  n-spectra  relations  via  minimum  cycle  basis 

2.1  Choose  decoupling  approximation  according  to  SNR 
regime  (c.f.  Proposition  3.3  or  3.4) 

2.2  Assign  baseline  edge  weights  accordingly 

2.3  Enumerate  the  Horton  cycles  of  the  interferometric 
graph 

2.4  Execute  greedy  algorithm  on  set  of  Horton  cycles  to 
find  minimum- variance  set  of  m  -  (Nap  -  1)  n-spectra 

3.  Average  the  selected  n-spectra 

4.  Unwrap  the  generalized  closure  phases  corresponding 
to  these  n-spectra  (c.f.  Section  3.1) 

5.  Solve  generalized  least-squares  problem  with  un¬ 
wrapped  generalized  closures 

6.  De-convolve  the  estimated  complex  visibilities  to  recon¬ 
struct  image  using  the  regularization  in  Equation  (65)  or 
other  regularization  technique 


estimation  in  the  absence  of  atmospheric  turbulence.  We  be¬ 
gin  by  defining  our  fringe  measurement  model  in  the  stan¬ 
dard  way  (see,  e.g.  Gordon,  J.  A.  &  Buscher,  D.  F.  (2012)): 


inary  parts  of  the  distinct  complex  visibilities  to  those  of  the 
generated  fringes.  The  FIM  then  becomes: 

K„d  =  (RrQrA^QR)_1  (69) 

Defining  the  function  gg  =  arctan  icv,d  ,  the  CRLB  for 

xcv,d 

the  phase  covariance  matrix  Cg  can  be  expressed  as: 

1  3ge(zCVjd)  -  <5gg(zcvj) 

Ce  Nf  dz  a  Izc'’-d  dz  a  ~  ° 

iy  frames  U£jcv,a  utjcv,a 

where  v-d)  is  the  Jacobian  matrix  of  the  distinct 

£*ZCV,£/ 

Fourier  phases  with  respect  to  the  distinct  complex- visibility 
phasors. 

Since  the  left-hand-side  of  Inequality  (70)  is  positive 
semi-definite  (PSD),  and  the  diagonal  entries  of  a  PSD  ma¬ 
trix  must  be  non-negative,  we  arrive  at  the  following  useful 
bound  for  the  variance  of  the  estimated  phases: 


var(di)  > 


1 

Nframes 


dgfl(zcv,<f)  i  dgQ{zCV'd)T 

dz cv,d  cx,d  dzCV(j 


ii 


(71) 


P  =  Q 


*cv 

ycv 


+  7/1 


(66) 


where  p  is  the  vector  of  pixel  counts,  xcv,  and  ycv  are  re¬ 
spectively  the  rea,l  and  imaginary  components  of  the 

complex-visibility  components,  respectively,  Q  is  the  matrix 
mapping  the  real  and  imaginary  parts  of  the  complex  visibil¬ 
ities  to  pixels  on  the  focal  plane  (often  called  the  visibility- 
to-pixel  matrix  V2PM  in  the  literature),  1  is  the  all-ones 
vector,  and  //  is  a  constant  scalar  whose  value  is  propor¬ 
tional  to  the  overall  brightness  of  the  target.  Note  that  rows 
of  Q  are  the  sinusoidal  functions  associated  with  each  fringe 
generated  by  the  beam  combiner. 

Leveraging  analysis  by  Zmuidzinas  (2003),  we  can  com¬ 
pute  the  Fisher-Information  Matrix  (FIM)  for  interferomet¬ 
ric  estimation  with  a  Fizeau  beam-combiner  as: 


Ixcv,ycv  =  (QrA^Q)"‘  (67) 

where  A fp  is  a  diagonal  matrix  whose  diagonal  entries  are 
the  expected  values  of  the  photon  counts  at  each  detector  in 
the  array. 

The  CRLB  for  the  covariance  of  these  vector  parameters 
is  then  given  by: 


>  0 


(68) 


where  the  notation  >  in  this  context  means  that  the  matrix 
difference  on  the  left-hand  side  is  positive-semidefinite. 

The  above  bound  applies  to  single-frame  fringe  phasor 
estimation  for  a  non-redundant  array.  Since  we  are  instead 
evaluating  phase  estimation  schemes  for  a  redundant  array 
given  multiple  frames  of  data,  we  must  incorporate  these 
transformations  into  the  bound  (Kay  1993).  Namely,  for  a 
redundant  array,  the  parameter  vector  must  be  shortened 
to  the  d  distinct  complex-visibility  phasors  xcv,d  +  (l/’)y cv,d- 
Let  us  define  the  vector  zcv  g  as  the  concatenation  of  xcvj 
and  yCVjj.  The  fringe  model  is  then  given  by:  p  =  QRzcv  j, 
where  R  is  a  l{N°p^  x2 d  matrix  mapping  the  real  and  imag- 


4.2  Simulation 

In  this  section,  we  provide  Algorithm  1  simulation  results 
for  both  interferometric  architectures.  For  these  simula¬ 
tions  we  used  an  RSC  pattern  of  the  Y-pattern  type  as 
shown  in  Figure  5.  The  corresponding  UV-sampling  for 
this  pattern  is  displayed  in  Figure  6.  To  prevent  alias¬ 
ing  on  the  simulated  focal  plane,  this  aperture  pattern 
was  mapped  onto  the  Golay  non-redundant  pattern  (Go- 
lay  1970)  shown  in  Figure  7  for  fringe  generation;  the 
model  emulates  the  redundant-to-non-redundant  pupil  re¬ 
mapping  technique  developed  by  Perrin  et  al.  (2006)  and 
described  in  the  Introduction.  Our  simulation  assumed 
Poisson-distributed  shot-noise  and  idealized  detectors  with 
zero  read  noise,  and  so  is  representative  of  a  shot-noise- 
dominated  scenario.  The  target  selected  was  NASA’s  Cloud- 
Aerosol  Lidar  and  Infrared  Pathfinder  Satellite  Observa¬ 
tions  (CALIPSO)  satellite  for  which  the  truth  image  (Hill 
2008)  is  shown  in  Figure  8,  as  well  as  a  reference  im¬ 
age  of  the  same  at  the  resolution  attainable  by  the  pat¬ 
tern.  Two  flux  levels  (2000  photoelectrons/aperture/frame, 
and  500  photoelectrons/aperture/frame)  were  considered  for 
these  simulations.  Comparative  error  analysis  using  Equa¬ 
tion  (57)  showed  a  lower  predicted  RMS  phase  error  for  the 
decoupling  approximation  in  Proposition  3.3  than  that 
in  Proposition  3.4,  and  hence  the  former  was  chosen  in 
Algorithm  1.  Bispectra  and  n-spectra  were  integrated  for 
50, 000  frames,  which  corresponds  approximately  to  an  8- 
minute  observation  time  if  we  assume  a  typical  frame  du¬ 
ration  of  10  msec.  An  implementation  of  the  Bellman-Ford- 
Moore  shortest-path  algorithm  (O’Connor  2012)  was  used 
to  generate  a  minimum  cycle  basis  as  per  Algorithm  1. 

To  show  the  potential  impact  of  generalizing  the  closure 
phase  notion,  Algorithm  1  was  compared  with  an  analogous 
algorithm  which  instead  used  the  minimum-variance  set  of 
^N“p  ^  independent  standard  (i.e.  three-baseline)  closure 
phases.  This  minimum  set  was  generated  via  application 
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Figure  5.  RSC  aperture  pattern  used  in  simulation 
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Figure  6.  UV-sampling  for  RSC  pattern 
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igure  7. 


Golay  non-redundant  beam-combiner  pattern 
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of  the  greedy  algorithm  described  above  to  the  set  of  all 

m- 

Numerical  results  for  the  pairwise  and  Fizeau  architec¬ 
tures  in  the  higher  flux  scenario  are  shown  in  Figures  9  and 
10,  respectively.  In  these  plots,  the  RMS  Fourier-phase  er¬ 
rors  for  10  independent  simulation  trials  are  plotted  as  a 
function  of  their  corresponding  visibilities  for  the  basic  so¬ 
lution.  Median-filtered  versions  of  the  predicted  standard 
errors  based  on  Equations  (57)  and  (61)  are  also  shown  for 
both  basic  and  pseudo-inverse  solutions,  respectively.  6 

As  expected,  the  Fizeau  architecture  proves  to  be  more 
sensitive  than  the  Pairwise  architecture,  which  corroborates 
analysis  by  Zmuidzinas  (2003).  Figures  11  and  12  respec¬ 
tively  show  Fizeau  performance  for  the  higher  flux  scenario 
at  a  shorter  integration  time  (1000  frames),  and  for  the  lower 
flux  scenario  at  the  long  integration  time  (50,000  frames). 
Note  that  the  use  of  independent  generalized  closure  phases 
according  to  Algorithm  1  outperforms  the  analogous  scheme 
in  which  only  standard  closures  are  used;  the  root-mean- 
squared  (RMS)  errors  are  noticeably  lower  for  the  former. 

It  is  interesting  to  examine  performance  relative  to  the 
CRLB  described  in  Section  5.1.  Recall  that  our  CRLB  is  an 
atmosphere-oracle  bound  and  hence  it  should  bound  phase 
estimation  accuracy  for  the  case  of  a  stationary  fringe  aver¬ 
aged  over  many  frames.  The  variance  of  such  an  estimator 
can  be  easily  calculated  by  leveraging  the  analysis  in  Sec¬ 
tion  2.  Namely  the  pseudovariance  of  the  phasor  is  given  by 
applying  Equation  (37)  after  averaging,  so  that  SNR  of  the 
averaged  fringe  phasor  z  is  given  by: 


1 

SNRi 


Nc 


ap 


I  T I Z 1 1  2N frcimestly 


(72) 


We  compute  the  predicted  phase  RMS  for  the 
atmosphere-oracle  by  taking  the  square-root  of  this  expres¬ 
sion,  and  add  it  to  the  plots.  From  the  plots,  we  see  that 
CRLB  is  virtually  identical  to  this  expression,  except  for 
isolated  dips  in  the  CRLB.  These  dips  correspond  to  the  re¬ 
dundant  baselines,  which  as  expected,  admit  increased  phase 
sensitivity  due  to  their  multiplicity. 

In  the  higher-flux  scenario  ( n  =  2000),  our  algorithm’s 
phase-estimation  performance  is  within  a  factor  of  2  from  the 
CRLB  for  the  vast  majority  of  sampled  visibilities.  As  the 
flux  drops  to  n  =  500,  the  performance  begins  to  diverge  from 
the  CRLB  as  the  closure  phase  variance  leaves  the  regime  de¬ 
scribed  by  Equation  (22)  and  enters  the  regime  described  by 
Equation  (21).  That  is,  the  variance  is  no  longer  accurately 
modeled  simply  by  a  linear  combination  of  the  individual 
phase  variances  of  the  form  in  Equation  (72);  rather  it  be¬ 
comes  inversely-proportional  to  the  product  of  the  squared- 
visibilities.  The  rapid  decrease  in  fidelity  of  the  closure  phase 
at  low-SNR  is  well-known  in  interferometry  (Kulkarni  et  al. 
1991). 

Sample  image  reconstructions  are  shown  in  Figure  13. 


6  The  actual  predictions  exhibit  oscillations,  which  reflect  the 
particular  structure  of  our  mapping  between  Fourier  phases 
and  GC’s.  For  clarity  we  have  created  smoothed  versions  of 
the  predicted  curves  by  creating  bins  with  edges  at  visibilities 
{5,  10,  20,40,  80,  160,  320 }  *  10  3 ,  and  displaying  the  median  values 
in  each  bin. 
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Figure  8.  Truth  image  for  simulation:  the  CALIPSO  satellite  (left),  and  Truth  image  at  the  resolution  of  the  interferometric  pattern 
(right) 


The  reconstructions  derived  from  generalized  closures  show 
greater  fidelity  to  the  true  image  than  those  derived  from 
standard  closures,  and  thereby  corroborate  the  numerical 
analysis  presented  in  the  plots.  To  quantify  relative  im¬ 
age  quality,  we  have  used  the  well-known  Struc¬ 
tural  Similarity  Metric  (SSIM)  (Wang  et  al.  2004). 
The  SSIM  computation  decomposes  the  image  into 
blocks  and  calculates  similarity  among  correspond¬ 
ing  blocks  in  the  test  image  and  the  reference  image. 
It  is  customary  to  report  the  mean  of  these  block- 
wise  similarity  scores,  i.e.  the  Mean  Structural  Simi¬ 
larity  Metric  (MSSIM).  The  value  of  MSSIM  ranges 
from  -1,  indicating  maximum  perceived  structural 
dissimilarity,  and  1,  indicating  perfect  structural 
similarity  (i.e.  identical  images).  The  MSSIM  score 
for  each  image  reconstruction  is  reported  in  its  cap¬ 
tion. 

It  is  noteworthy  that  while  our  algorithm  attempts  to 
find  a  minimum- variance  set  of  n-spectra,  this  set  is  not  nec¬ 
essarily  the  optimum  set  for  phase  estimation;  our  algorithm 
does  not  account  for  the  conditioning  of  the  matrix  C 0c  re¬ 
sulting  from  a  particular  choice  of  nspectra.  In  fact  there 
may  be  some  cases  in  which  the  use  of  traditional  closures 
results  in  lower  phase-estimation  error  due  to  better  condi¬ 
tioning  in  the  corresponding  Coc.  In  practice,  the  decision 
to  use  traditional  or  generalized  closures  can  be  resolved 
by  examining  the  phase  errors  predicted  by  Equation  (57), 
for  basic  solutions,  or  Equation  (61)  for  pseudo-inverse  so¬ 
lutions. 


4.3  Generalized  Closures  in  Non-Linear  Least 
Squares  Approaches 

This  paper  has  proposed  both  novel  reconstruction  method¬ 
ology  as  well  as  new  interferometric  observables.  A  natural 
question  arises  as  to  whether  the  advantages  of  the  latter  are 
retained  with  other,  existing  methodology.  In  particular,  we 
consider  the  class  of  techniques  which  solve  the  non-linear 
least  squares  inference  problem  involving  normalized  bis¬ 


pectrum  phasors  (Gorham  et  al.  1989),  (Negrete- Regagnon 
1996).  Such  approaches  minimize  an  objective  of  the  form: 


|L jfii.ct  —  ej(6n+8i2+9n)\\ 


„,2 


(73) 


where  Nc  is  the  number  of  closures,  is  the  z'-th  closure 
phase,  the  Gy,  8/2,  and  8/j  are  the  Fourier  phases  associated 
with  this  closure,  and  w/  is  a  weighting  factor  proportional 
to  the  estimated  variance  of  the  closure  measurement.  Ob¬ 
jective  T  is  clearly  a  continuous  function  and  hence  local 
minima  can  be  reliably  found  with  gradient-descent  tech¬ 
niques. 

While  such  algorithms  bring  the  complexities  involved 
in  non-linear  optimization  (e.g.  possible  stagnation  in  lo¬ 
cal  minima  and/or  slow  convergence),  they  also  feature  the 
ability  to  utilize  any  number  of  closure  relations.  Algorithm 
1  restricts  attention  to  closure  relations  which  form  cycle 
bases,  thereby  guaranteeing  well-posed  Fourier  phase  recov¬ 
ery  while  keeping  the  size  of  the  associated  CVP-unwrapping 
problem  tractable.  However,  as  noted  by  Kulkarni  et  al. 
(1991),  fitting  to  all  closures  will  in  principle  improve 

estimation  performance  when  the  per-frame  flux  is  low,  due 
to  the  fact  that  closure  phases  de-correlate  as  the  flux  de¬ 
creases  7.  Specifically,  as  computed  by  Kulkarni  et  al.  (1991) 
for  traditional  closures  in  the  pairwise  case,  the  correlation 
coefficient  y  of  two  closures  sharing  a  common  baseline  ap¬ 
proaches  ^  in  the  high-SNR  regime  (i.e.  for  hy~  »  1)  as  one 
would  expect  8.  On  the  other  hand,  the  correlation  coeffi¬ 
cient  in  the  low-SNR  regime  can  be  accurately  modeled  as 
y  ~  4 y4n  • 

To  assess  the  phase-estimation  performance  provided 
by  generalized  closures  relative  to  that  using  the  complete 


1  The  interested  reader  is  directed  to  Buscher  (2015)  for  a  de¬ 
tailed  discussion  of  this  issue 

8  Here  Kulkarni  et  al.  (1991)  assumes  a  common  visibility  y 
among  all  baselines  in  the  closure 
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Figure  9.  Pairwise  Phase  Recovery  Results  for  Flux  n  =  2000  pe/ap/frame,  50,000  frames 
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Figure  10.  Fizeau  Phase  Recovery  Results  for  Flux  n  =  2000  pe/ap/frame,  50,000  frames 
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Figure  11.  Fizeau  Phase  Recovery  Results  for  Flux  n  =  2000  pe/ap/frame,  1000  frames 
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Figure  12.  Fizeau  Phase  Recovery  Results  for  Flux  n  =  500  pe/ap/frame,  50,000  frames 
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Figure  13.  Fizeau  Image  Reconstruction  Results.  The  reconstructions  in  the  left  column  used  traditional,  three-baseline  observ¬ 
ables,  whereas  those  in  the  right  column  used  generalized  closures  selected  according  to  Algorithm  1.  (top  row)  n  =  2000  photoelec¬ 
trons/aperture/frame  (pe/ap/frame),  50,  000  frames,  (middle  row)  n  =  2000  pe/ap/frame,  2000  frames,  (bottom  row)  n  =  500  pe/ap/frame, 
50,  000  frames. 


MSSIM  0.8466  MSSIM  0.8646 


MSSIM  =  0.5940  MSSIM  =  0.6841 


MSSIM  =  0.6623  MSSIM  =  0.7166 
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set  of  traditional  closures,  we  consider  the  algorithm  devel¬ 
oped  first  by  Gorham  et  al.  (1989)  which  works  the  normal¬ 
ized  bispectrum  instead  of  the  closure  phase.  We  used  the 
MATLAB®  non-linear  least-squares  solver  (NLS)  known  as 
Isqnonhn  (MATLAB  2014)  to  minimize  the  objective  db  in 
Equation  (73).  This  solver  obtains  a  quadratic  approxima¬ 
tion  of  the  objective  at  each  step  and  moves  towards  the  min¬ 
imum  of  this  approximation,  at  which  point  another  approx¬ 
imation  is  computed  and  the  process  repeats.  The  weights 
in  db  were  set  as  in  Gorham  et  al.  (1989): 


where  ||g,jl  is  the  magnitude  of  the  averaged  unbiased  nspec- 
trum,  and  a-gi  is  the  empirical  standard  deviation  of  its 
quadrature  components. 

Results  of  applying  this  solver  to  the  lower-SNR  sce¬ 
nario  of  the  previous  section  are  given  in  Figure  14.  A  single 
randomly-chosen  initialization  point  was  used  in  all  cases. 
The  top  row  shows  reconstructions  and  associated  conver¬ 
gence  times  for  the  NLS  algorithm  using  traditional  closures 
(left)  and  generalized  closures  selected  from  a  minimum  cy¬ 
cle  basis  (middle),  respectively.  The  bottom  row  gives  anal¬ 
ogous  results  for  the  case  in  which  all  (  traditional  clo¬ 
sures  are  used.  Three  snapshots  at  iteration  counts  20,  40, 
and  80  are  shown.  Iterations  after  80  resulted  in  negligible 
change  to  the  image  quality.  These  results  suggest  that  gen¬ 
eralized  closures  can  provide  at  least  the  sensitivity  of  tra¬ 
ditional  closures.  At  the  same  time,  they  benefit  from  the 
increased  convergence  speed  afforded  by  a  problem  size  pro¬ 
portional  to  the  size  of  cycle  basis  (i.e.  <x  N^p)  as  opposed  to 
the  size  of  the  set  of  all  closures  (i.e.  oc  N^p).  These  results 
hence  show  that  generalized  closures  can  serve  as  more  effi¬ 
cient  sources  of  phase  information  than  traditional  closures. 

In  Figure  15,  we  plot  phase  error  as  a  function  of  visi¬ 
bility  for  the  NLS  approach  using  a  cycle  basis  of  traditional 
closures  (blue),  a  cycle  basis  of  generalized  closures  (green), 
and  using  all  (  traditional  closures  (red).  The  median 
absolute  error  for  20  trials  is  shown  in  the  left  panel.  The 
same  median-filter  that  was  applied  to  the  predictions  in 
Figures  9-12  was  then  applied  to  this  raw  data  to  produce 
the  results  in  the  right  panel  of  the  Figure.  It  is  clear  that 
the  sensitivity  advantage  of  generalizing  the  closure  relations 
are  retained  in  the  NLS  approach  for  this  example. 

It  is  noteworthy  that  the  selection  of  a  minimum- 
variance  set  of  generalized  closures  involves  computational 
overhead  in  the  execution  of  the  greedy  algorithm.  How¬ 
ever,  it  is  anticipated  that  in  many  cases  we  may  be  able 
to  obtain  a  reliable  surrogate  for  this  set  by  instead  form¬ 
ing  a  fundamental  cycle  basis  (Liebchen  &  Rizzi  2005)  as¬ 
sociated  with  minimum-variance  spanning  tree  of  the  inter¬ 
ferometric  graph.  Constructing  this  alternate  basis  obviates 
the  need  for  execution  of  the  greedy  algorithm;  it  merely 
requires  the  execution  of  a  Minimum-weight  Spanning-tree 
algorithm  (e.g.  Prim-Jarnik  Algorithm  (Prim  1957))  with 
each  edge  weight  set  to  reciprocal  of  the  squared  visibility 
associated  with  baseline  9.  The  elements  of  the  correspond- 

9  For  our  simulation,  we  used  a  MATLAB®  implementation  of 
Prim’s  algorithm  which  is  available  online  (Greenbaum  2007) 


ing  minimum- variance  spanning  tree  (MVST)  cycle  basis  are 
then  found  by  closing  the  edges  of  this  spanning  tree,  which 
is  shown  in  Figure  16.  Indeed,  many  of  the  Horton  cycles 
comprising  the  minimum  cycle  basis  computed  for  this  ex¬ 
ample  pattern  and  scenario  are  also  members  of  the  MVST 
basis.  Not  surprisingly,  the  reconstruction  with  this  basis 
approaches  the  quality  of  that  for  the  minimum-cycle-basis 
set  (c.f.  middle  and  bottom  images  in  left  column  of  Figure 
14).  Moreover  the  median  absolute  phase  error  for  this  ap¬ 
proach,  which  is  shown  in  cyan  in  the  right  panel  of  Figure 
15,  is  very  similar  to  that  observed  with  the  Minimum  Cycle 
Basis. 

While  the  analysis  of  Kulkarni  et  al.  (1991)  revealed 
the  advantage  of  supplementing  a  cycle  basis  with  addi¬ 
tional  traditional  closures  in  the  pairwise  case,  this  advan¬ 
tage  should  in  principle  extend  to  the  case  of  generalized 
closure  phases  in  the  Fizeau  architecture.  Our  CRLB  com¬ 
parison  above  (i.e.  Figures  10  and  11)  confirms  that  there  is 
limited  scope  for  improvement  in  the  regime  jp—  »  60.  In 
the  regime  explored  in  the  lower-SNR  scenario  (see  Figure 
12),  however,  there  is  certainly  a  CRLB-gap  which  may  be 
closed  with  appropriate  selection  of  supplementary  gener¬ 
alized  closures.  We  leave  the  development  of  algorithms  to 
perform  such  a  selection  efficiently  for  future  work. 


5  CONCLUSION 

In  this  paper,  we  have  developed  a  novel  method  for  interfer¬ 
ometric  imaging  which  employs  RSC  techniques  and  a  gener¬ 
alized  notion  of  bispectrum  observable  (the  n-spectrum) .  We 
have  established  that  phase  estimation  from  these  observ¬ 
ables  is  well-posed  for  valid  RSC  arrays.  We  have  also  pro¬ 
vided  a  fast  algorithm  for  selection  of  a  minimum-variance 
set  of  these  observables,  which  is  based  on  the  concept  of 
minimum  cycle  basis  in  graph  theory.  Then,  leveraging  the 
lattice-theory  problem  formulation  first  proposed  by  Lannes 

6  Anterrieu  (1999)  for  unwrapping  of  the  closure  phases, 
as  well  as  techniques  from  sparse  recovery  for  image  re¬ 
construction  in  the  presence  of  Fourier  undersampling,  we 
have  proposed  a  new  algorithm  for  image  reconstruction 
in  optical  interferometry.  Previous  RSC-based  algorith¬ 
mic  frameworks  have  used  the  measured  baseline 
phase  as  the  primary  observable  (Lannes  &  Anter¬ 
rieu  1999)  and/or  rely  on  computation  of  the  Smith 
Normal  Form  (Lannes  2003).  Our  framework  in¬ 
stead  operates  on  those  atmosphere-invariants  avail¬ 
able  in  a  typical  long-baseline  interferometric  imag¬ 
ing  scenario,  and  leverages  standard  least-squares 
techniques  tailored  to  the  statistics  of  these  ob¬ 
servables.  We  have  confirmed  that  the  performance  of  the 
phase-estimation  part  of  our  algorithm  can  be  easily  quan¬ 
tified  from  first  principles  using  standard  linear  estimation 
theory.  Finally,  we  have  conducted  simulations  at  differ¬ 
ent  photon- flux  levels  to  corroborate  this  analysis,  and  to 
show  the  potential  advantage  of  performing  inference  on  the 
n-spectrum  observable  relative  to  its  classical  counterpart 
-  the  bi-spectrum.  Large-scale  optical  systems  such  as 
the  Planet  Formation  Imager  are  currently  being  de¬ 
signed  to  support  ambitious  missions  including  ex¬ 
oplanet  observation.  It  is  our  hope  that  both  the  the¬ 
oretical  framework  employed  in  this  paper,  as  well  as  the 
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Figure  14.  (top  row)  Reconstructed  images  and  convergence  times  for  NLS  algorithm  with  traditional-closure  basis  (left),  generalized- 
closures  from  minimum  cycle  basis  (middle),  and  generalized-closures  from  minimum- variance  spanning  tree  (bottom  row)  Reconstructed 
images  and  elapsed  running  times  with  all  p )  traditional  closures  at  iteration  20  (left),  iteration  40  (middle),  and  iteration  80  (right) 


Cycle  Basis  (Traditional)  All  Traditional  Closures 

MSSIM  =  0.625,  Convergence  in  5  sec  MSSIM  =  0.2175,  20  iterations,  73  sec  elapsed 


Cycle  Basis  (Generalized,  Min  Cycle  Basis) 
MSSIM  =  0.686,  Convergence  in  4  sec 


All  Traditional  Closures 
MSSIM  =  0.6195,  40  iterations,  138  sec  elapsed 


Cycle  Basis  (Generalized,  Min.  Spanning  Tree) 
MSSIM  =  0.750,  Convergence  in  3  sec 


All  Traditional  Closures 
MSSIM  =  0.6133,  80  iterations,  258  sec  elapsed 
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Figure  15.  Raw  (left)  and  median-filtered  (right)  Fourier  phase  error  using  the  NLS  approach. 


Figure  16.  The  minimum- variance  spanning  tree 


practical  algorithm  itself,  will  prove  useful  to  future  users  in 
designing  such  systems  and  processing  their  data. 
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APPENDIX  A:  FIZEAU  VARIANCE 
APPROXIMATIONS 

Al  Variance  Decomposition 

In  this  section  we  justify  our  approximation  for  the  variance 
of  the  n-spectrum  for  the  Fizeau  architecture.  Our  goal  here 
is  not  to  compute  all  of  the  numerous  terms  in  the  variance, 
but  rather  to  present  the  mathematical  intuition  behind  the 
approximation  we  have  used.  The  general  formula  for  the 
variance  of  an  n-spectra  G  is: 


var[G] =(nz>'4)-  (nuHn  4)  (ai) 

'i=1  I  '1=1  I  'f=i  I 

For  simplicity,  we  analyze  the  standard  bispectrum  case 
(o  =  3).  The  analysis  extends  naturally  to  the  n-spectrum. 
Note  that  by  extension  of  Equation  (5)  the  first  term  can  be 
written  as: 


1” Mr  X  {q{pa)q{pb)q(pc)q{pd}q(pe)q{pf)) 

1=1  '  Pa,Pb,Pc,Pd,Pe,Pf 

x  ei<bi(Pd-Pa)  giO>2(pe-Pb)  gUOiipf-pa)  (A2) 

Similarly,  the  second  term  in  Equation  (Al)  can  be  writ¬ 
ten  as: 

([>)([>  ?)  = 

h‘=  1  '  *i=l  ' 

^  < q(Pa)q(pb)q(Pc))(q(Pd)q(Pe)q(pf )> 

Pa,Pb,Pc,Pd,Pe,Pf 

x  giMl(pd-Pa) gi<P2(Pe-Pb) gi<Pl(Pf~Pc)  (A2>) 

Proceeding  analogously  to  Equation  (7),  we  utilize  the 
general  formula  for  the  moments  of  a  Poisson  distribution 
to  perform  the  following  decomposition  of  the  first  factor  in 
the  summand  in  Equation  (A2)  (Kulkarni  et  al.  1991): 
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(q(Pa)q(Pb)q(Pc)q(Pd)q(Pe)q(Pf))  =  (A4) 

+(q(pa))(q(pb))(q(pc))(q(pd))(q(pe))(q(pf))  (A5) 

+6p=pa=pj(q(p))(q(pb))(q(Pc))(q(Pe))(q(Pf))  (A6) 

+Sp=Pb=Pe  (q(p))(q(Pa))(q(Pc))(q(Pd))(q(Pf))  (A7) 

+6p=Pc=pf(q(p))(q(pa))(q(pb))(q(Pd))(q(Pe))  (AS) 

+3p=pa=Pd(q(p))6p'=pb=Pe(q(p'))(q(pc))(q(pf))  (A9) 

+Sp=Pa=Pd{q(p))6p'=Pc=Pf{q(p')){q(pb)){q{pe )>  (A10) 

+Sp=Pb=Pe(q(p))6p'=Pc=Pf(q(p'))(q(pa))(q(pc ))  (All) 

+Sp=Pa=Pd(q(p))Sp'=Pb=pe(q(p  ))Sp"=Pc=Pf(q(p  ))  (A12) 

+8p=Pa=Pe(q(p))(q(pb))(q(pc))(q(pd))(q(Pf))  (A13) 

+8p=Pa=Pf(q(p))(q(pb))(q(pc))(q(Pd))(q(Pe))  (A14) 

+6p=pb=Pd(q(p))(q(pa))(q(Pc))(q(Pe))(q(Pf))  (A15) 

+8p=Pb=Pf(q(p))(q(pa))(q(Pc))(q(Pd))(q(Pe))  (A16) 

+8p=Pc=Pd{q(p)){q(pa)){q{pb)){q(Pe)){q(Pf))  (A17) 

+8p=Pc=Pe(q(p)){q(pa))(q(Pb))(q(Pd))(q(Pf))  (A18) 

+8p=Pa=Pb(q(p)){q(Pc))(q(Pd))(q(Pe))(q(Pf))  (A19) 

+8p=Pa=Pc(q(p))(q(pb))(q(Pd))(q(Pe))(q(Pf))  (A20) 

+8p=Ph=Pc(q(p))(q(pa))(q(pd))(q(pe))(q(Pf))  (A21) 

+8p=Pd=Pe(q(p))(q(pa))(q(Pb))(q(Pc))(q(Pf))  (A22) 

+Sp=Pd=Pf(q(p))(q(pa))(q(pb))(q(Pc))(q(Pe))  (A23) 

+8p=Pe=Pf(q(p))(q(pa))(q(ph))(q(Pc))(q(Pd))  (A24) 

+8p=Pa=Pb=Pc(q{p))(q(pd))(q(Pe))(q{pf ))  (A25) 

...  +  other  third-order  terms  (A26) 

+8p=Pa=Pb=Pc=Pd(q(p))(q(pe))(q(pf))  (A27) 

...  +  other  fourth-order  terms  (A28) 

...etc.  (A29) 


2  <«(?)><« (Pb))eioJ2Pb  {q{Pc))ei(wi+on)^ 

P=Pa=Pd,Pb,Pc,Pe,Pf 

-  Z  ap  Z  eiW2Pfa  Z  A^  e''(w,+"2)p<; 

p.—Pa  ~Pd  Pb  Pc 

Z  APeeiW2Pc  ZAP/e”'<Wl+"2)P/ 

Pe  Pf 

(A30) 

The  magnitude  of  this  expression  is  (Napn-’y^y^).  In 
the  case  where  all  visibilities  are  equal  to  y.  this  expression 
becomes  Napn5y 4.  Terms  (A7)-(A8)  follow  this  same  form. 
Now  consider  term  (A9): 

Z  {q(p)){q(p')){q{pe))e-ia)2P‘ 

P=Pa  =Pd,P'=Pb  =Pe,Pe,Pf 

x  {q(pf))e-i(0il+COl)Pf 

=  Z  AP  Z  ApZA  Pce^+^)Pf 

P-=Pa=Pd  P'-=Pb=Pe  Pc 

x  Yj  APfe~ii0Jl+a>2)pf  (A31) 
Pf 

The  magnitude  of  this  expression  is  N^pn^y^,  which, 
under  the  assumption  of  equal  visibilities,  becomes  N^pn4y2. 
Terms  (AlO)-(All)  follow  this  same  form. 

Case  2:  mixed  non-conjugate  pairs 

We  now  consider  the  terms  in  which  a  pair  includes  two 
distinct  baselines  (e.g.  (A13)-(A18)).  The  common  pixels  in 
this  case  are  such  that  one  of  the  corresponding  baselines  is 
conjugated  and  the  other  is  not.  Taking,  for  example,  the 
case  (A13)  in  which  p  :=  pa  =  pe,  we  have: 


The  total  number  of  partitions  (203)  in  the  summation 
above  is  given  by  the  6th  Bell  number.  Below  we  categorize 
the  terms  by  their  order,  which  in  this  case  we  define  as  the 
largest  number  of  common  pixels  in  each  partition.  We  will 
see  that  we  can  approximate  the  pseudo-variance  to  reason¬ 
able  accuracy  by  including  a  particular  subset  of  the  Order-2 
terms.  Note  that  all  terms  in  the  analogous  decomposition 
of  Equation  (A3)  are  also  in  Equation  (A2)  and  hence  cancel 
them  in  Equation  (Al). 


A2  Order-2  Terms 

These  terms  have  a  pair  of  common  pixels.  We  distinguish 
between  the  following  two  cases: 

Case  1:  conjugate  pairs  These  terms  are  given  in  (A6)- 
(A12).  The  common  pixels  belong  to  the  conjugate  pairs,  i.e. 
P  '■=  Pa  =  Pd,  P  '■=  Pb  =  Pe,  and  p  :=  pc  =  pf,  respectively, 
which  also  appear  in  the  pairwise  case.  Let  us  assume  that 
the  visibilities  associated  with  a>\,  a>2,  and  W3  are  all  roughly 
equal,  i.e.  y  =  y \  =  72  =  73-  Consider  the  sum  of  the  first  of 
these  terms  (A6)  over  the  focal  plane,  i.e.: 


^  (q(p))ei^-^p(q(pb))ei‘P2P» 

P=Pa  =Pe,Pb,Pc,Pd,Pf 

x<?(Pc)>ei(Wl^)Po<®(Pd)>e-^<9(P/)>e-<(Wl+Wl)^ 

=  Y  ^pei<w,-°J2)p  Z  Z  V  e'("'+“2)P' 


P  —P  a =Pe 


Pb 


xZAP/"iW  Zv  e-i(.toi+tJ2)pf  (A32) 


Pd 


Pf 


where  y0)]  is  the  visibility  of  the  fringe  of  spatial  fre¬ 
quency  a>  1  -  u>2 .  Note  that  the  first  factor  in  the  Equa¬ 
tion  above  clearly  vanishes  if  such  a  fringe  is  not  created 
by  the  beam  combiner.  Otherwise,  assuming  equal  strength 
of  the  visibilities,  the  magnitude  of  this  expression  is  given 
by  (n5 yWl-a>2y\y2y2) ,  or  n5y5  in  the  case  of  equal  visibilities. 
Comparing  with  any  of  the  Case  1  terms  and  recalling  that 
y  <  1 ,  we  see  that  the  latter  will  dominate  the  former  in  the 
regime  considered  in  this  paper  («  as  le3,  y  1,  Nap  =31). 

Case  3:  terms  in  decomposition  Equation  (A3) 

This  Case  includes  terms  (A5),  as  well  as  (A19)-(A24). 
The  common  pixels  are  such  that  both  or  neither  of  the 
corresponding  baselines  are  conjugated,  i.e.:  the  cases  pa  = 
Pb,  Pa  —  Pc ,  Pb  ~  Pc  Pd=Pe,Pd=  Pf,  Pe  =  Pf  -  These  terms 
are  canceled  by  their  counterparts  in  Equation  (A2).  □ 
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A3  Higher  Order  Terms 

An  example  of  a  third-order  term  would  be  the  case  of  p  :  = 
Pa  -  Pb  -  Pdi  whose  sum  is  given  by: 


P'=Pa  =Pb  =Pd,Pc,Pe,Pf 

(qiPeVe-^P'iqipfpe-^+^Pf 


(A33) 


=  ^  A p  ei0J2P  J]  A Pc  ei(LJr +"2  )Pc 

P-=Pa=Pb=Pd  Pc 

Ev  e-i(wi+w2)p/  (A34) 
Pe  Pf 

The  magnitude  of  this  expression  is  given  by  n^y^yf 
For  equal  visibilities,  we  have  «4y4.  Clearly  this  term  will 
be  dominated  by  the  Case  1  terms.  Other  high-order  terms 
exhibit  a  similar  sharp  attenuation  which  allows  them  to  be 
neglected  for  practical  purposes.  □ 

Given  the  relative  strengths  of  the  terms  shown  above, 
we  will  retain  Case  1  terms  (A6)-(A12),  which  after  sum¬ 
mation  and  factoring,  yield  Equation  (38).  Note  that  this 
amounts  to  the  following  approximation: 


Vfizeait(G)  '■ 


1=1 


f]<z;z,*>  -  f”[  > 


(A35) 


i=i 


Applying  these  same  approximations  leads  to  the  Fizeau 
covariance  expressions  in  Equations  (40)-(41). 


Corollary  B.2:  For  a  valid  RSC  array,  the  mapping  Coc  is 
injective  up  to  an  image  shift. 

Proof:  Proposition  B.l  showed  that  the  nullspace  is 
comprised  of  linear  combinations  of  vectors  Ar*  and  Ary. 
Note  that  Ar*  and  Ary  are  simply  scaled  versions  of  the  x- 
and  y-spatial  frequency  vectors  in  the  array,  respectively. 
Hence  adding  linear  combinations  of  these  vectors  to  a  par¬ 
ticular  RSC  phase  solution  merely  produces  phase  ramps  in 
the  Fourier  domain,  which  are  equivalent  to  translations  (or 
shifts)  in  the  image  domain.  In  other  words,  the  mapping 
C0c  is  invertible  up  to  an  unknown  image  shift.  □ 


APPENDIX  C:  PROOF  OF  PROPOSITION  3.2 

We  begin  the  Proof  with  the  following  Lemma: 

Lemma  C.l:  The  final  2  columns  of  Vo-  form  a  basis  for 
the  nullspace  of  Coc. 

Proof:  This  follows  from  the  fact  Coc  is  rank-deficient  by  2, 
and  standard  properties  of  the  right  singular  vectors  com¬ 
prising  Vo-  in  the  SVD.  (Bretscher  2001)  □ 

Note  that  the  error  resulting  from  application  of  the 
pseudo-inverse  C„c  to  the  unwrapped  vector  of  closures  will 
be  given  by: 


2nea  =  C+oc(2ne*h)  (Cl) 

Let  us  express  the  spatial  frequencies  measured  by  an 
array  as  two-element  vectors  of  the  form  (cox,ojy).  Let  X  be 
the  d  x  2  matrix  containing  these  spatial  frequencies.  Note 
then  that  the  phase-wrap  error  will  manifest  itself  merely 
as  an  image  shift  if  and  only  if  this  error  is  a  (modulo-2^) 
phase  ramp,  i.e.  there  exists  a  2-element  shift  vector  z  and 
an  integer  vector  k  which  together  satisfy: 


APPENDIX  B:  PROOF  OF  INJECTIVITY  OF 
THE  RSC  CLOSURE  MAPPING 

Proposition  B.l:  For  a  valid  RSC  array,  the  columns  Arf 
and  Ary  form  a  basis  for  the  two-dimensional  nullspace  of 

Coc 

Proof:  To  see  that  the  two  columns  mentioned  belong 
to  the  nullspace  of  C oc,  note  that  each  solution  set  to  the 
noiseless  version  of  Equation  (46)  above  remains  valid  after 
replacing  each  0,-y  with  BK.  =  Ojj  -  z  •  (r,-  -  r j). 

We  then  need  to  establish  that  these  two  vectors  span 
the  entire  nullspace  of  C oc  by  showing  that  this  nullspace 
is  two-dimensional.  Suppose  we  have  a  vector  w  which  is  in 
the  nullspace  Coc.  This  is  equivalent  to  the  either  of  the 
following  conditions:  M#w  =  0,  or  Mgw  s  kerf Cmc).  The 
former  condition  is  not  possible  since  the  columns  spanning 
the  subspace  K  are  linear ly-independent.  The  latter  condi¬ 
tion  is  equivalent  to  the  condition  that  Mjw  6  K  n  L.  It  is 
well-known  fact  that  for  any  two  subspaces  K  and  L,  we 
have: 

dimfK  n  L)  =  dim(K)  +  dim(L)  -  dimfK  +  L)  (Bl) 

We  established  in  Section  3.1  that  dimfK  +  L)  = 
rankf M)  =  d+N-3  for  a  valid  RSC  system.  Also  dimfK )  =  d, 
since  the  K  is  spanned  by  d  linearly-independent  columns. 
Finally  dim(L)  =  N  -  1.  Substituting  into  Equation  (Bl),  we 
see  that  dimfK  fll)  =  2.  □ 


2 necr  -  2nXz  =  2nk  (C2) 

Substituting  from  Equation  (Cl)  we  obtain: 

C+c(27re)  -  2 nXz  =  2 nk  (C3) 

Dividing  through  by  2 n  we  obtain  the  equation:  C^c  e*h  - 
Xz  =  k.  Note  that  each  element  of  C+c  can  be  expressed 
as  some  rational  number  Similarly  we  first  assume  X 
contains  rational  spatial  frequencies  with  greatest  common 
denominator  qx.  Then  we  can  multiply  through  by  the  least- 
common-multiple  (LCM)  of  the  {qp  and  qx  to  obtain  a  sys¬ 
tem  of  equations  whose  coefficients  are  guaranteed  to  be 
integer  (i.e.,  we  have  a  linear  Diophantine  system).  Let  this 
LCM  be  denoted  as  I.  Then  we  have,  after  rearranging  terms, 

IXz  =  /(C+ce*  -  k)  (C4) 

We  now  wish  to  determine  conditions  under  which  there 
exist  vectors  k  and  z  satisfying  this  overdetermined  Dio¬ 
phantine  system.  Applying  the  Smith  Normal  Form  decom¬ 
position  (c.f.  Theorem  3.1)  to  the  matrix  IX  this  time,  and 
noting  that  rankfX )  =  2,  we  have: 

DX  =  Ux(/X)Vx  (C5) 

where  Ux  and  Vx  are  unimodular  matrices  of  size  m  x  m  and 
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2x2,  respectively,  and  Dx  is  a  rectangular  diagonal  matrix 
whose  entries  are  zero  below  row  2. 

If  we  left-multiply  Equation  (C4)  by  Ux  on  both  sides, 
we  obtain: 

!UxXz  =  lUx(C+ce*h-k)  (C6) 

Using  Equation  (C5)  and  the  fact  that  Vx  is  a  unimod- 
ular  (and  hence  invertible)  matrix,  we  can  then  write: 

DXVx'z  =  /(UxC+ce*  -  Uxfc)  (C7) 

We  are  now  in  position  to  prove  the  main  result  of  this 
section,  which  is  preceded  by  the  following  Lemma: 
Lemma  C.2:  Given  wrap-invariance,  the  (column)  vector 
UXC+ce  has  integer  entries  below  row  2. 

Proof : 

Given  that  the  elementary  divisors  of  C 0c  are  all  1,  we 
know  there  exists  an  integer  vector  &q  such  that: 

e*  =  Cocfc0  (C8) 

Substituting  Equations  (C8)  and  (60)  into  Equation 
(Cl),  we  obtain: 

eo-  =  U^Uo-E^fco  (C9) 

Noting  that  U a-  is  orthogonal,  this  equation  can  be  sim¬ 
plified  to 

e£r  =  (\Za-  -  N)V£fc0  (CIO) 

where  N  is  a  matrix  of  the  same  size  as  V^,  which  is  zero 
except  for  the  last  two  columns.  These  last  two  columns  are 
identical  to  those  of  V^-,  and  hence  by  Lemma  4.1  comprise 
an  orthogonal  basis  for  the  nullspace  of  Coc.  Noting  the 
orthogonality  of  V^-,  this  can  be  further  simplified  to: 

eo-  =  fc0-NV^fc0  (Cll) 

Note  that  each  of  the  two  non-zero  column  vectors 
jtq),  t  6  1,2  in  N  can  be  expressed  as  a  linear  combina¬ 
tions  of  the  elements  of  the  canonical  basis  for  the  nullspace 
of  Coc  given  in  Section  3,  i.e. 

Vfr  =  +  C>2W2  (C12) 

Hence  we  have: 

Uxeo-  =  UxkQ  -  UxN^V^fco  (C13) 

Since  Ux  is  an  integer  matrix,  the  first  term  is  clearly 
integral.  Let  us  then  examine  the  second  term,  and  in  par¬ 
ticular,  the  product  UxN.  By  substitution  from  Equation 
(C12),  we  see  that  or  any  of  the  two  non-zero  columns  of 
N,  we  have: 

Uxi’/f  =  fliUxu’i  +  a2Uxul2  (C14) 

Now  note  from  Theorem  3.1  that  Ux  is  a  matrix  which 
annihilates  all  the  spatial  frequencies  in  the  matrix  X  below 


row  2.  But  recall  from  Section  3  that  these  spatial  frequen¬ 
cies  are  the  contents  of  the  two  columns  uq  and  uq  (up 
to  a  uniform  scaling  factor).  Therefore  the  column  vector 
in  Equation  (C14)  is  zero  below  row  2.  This  means  in  turn 
that  the  second  term  in  Equation  (C13)  is  zero  below  row 
2,  and  hence  that  UXe<r  is  integral  below  row  2  (since  the 
first  term  is  integral).  □ 

To  utilize  Lemma  C.2,  we  first  re-arrange  the  Equation 
(C7)  above  so  that  it  reads: 

}dxVx'z  -  \JxC+oce*h  =  -Uxfe  (C15) 

Let  v  =  jDxV^'z  -  UxC+ce*.  Note  that  since  Dx  is 
zero  below  row  2,  the  entries  of  v  below  row  2  will  be  equal 
to  those  of  (-UxC+ce*),  which  are  integers  by  Lemma  C.2. 
Now  consider  the  first  and  second  entries  of  v.  Let  /  be 
the  vector  containing  the  fractional  parts  of  the  first  two 
elements  of  vector  uxc^c  e*h,  and  let  A  be  the  invertible 
matrix  consisting  of  the  first  two  rows  of  j-DxV^1.  Without 
loss  of  generality,  choose  z*  -  A ~lf  so  that  the  fractional 
part  /  is  annihilated,  leaving  only  integer  elements  in  the 
first  two  entries  of  v.  Hence  we  now  have: 

w  =  -Uxfc  (C16) 

with  v  ensured  to  contain  only  integer  elements.  Since 
Ux  is  unimodular,  the  vector  k*  =  will  be  integral. 

We  have  thus  found  a  pair  (z*,  k*)  with  integer  k*  which 
satisfies  the  Equation  (C7).  Since  Equation  (C7)  is  related 
to  Equation  (C4)  via  a  unimodular  (and  hence  invertible) 
mapping  Ux,  invariance  is  hence  proven.  □ 
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